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SuMMARY: By a slight modification to den Hartog’s classical theory of 
harmonically forced vibration off resonance it is shown that, with such 
distorted forcing wave forms as might be encountered in the stalling or 
choking flutter of an aerofoil in cascade, the influence of damping on the 
amplitude of motion can still be assessed with considerable accuracy by the 
assumption of an equivalent harmonic excitation. The result is valid for any 
amount of Coulomb damping at which motion is significant in the presence 
of any “velocity” damping coefficient up to that defined by a logarithmic 
decrement of about 1-3, the formal conclusion of the argument being stated 
in Section 4. 


1. Introduction 


Among the stress components experienced by the blading of an axial com- 
pressor one of the most important is the periodic bending stress produced by such 
oscillations as are due partly to mechanical forcing through the disc or casing, 
partly to aerodynamic forcing by periodically incident wakes, and in the worst 
practical instances to self-induced aeroelastic oscillation of the blades involved. 
Such self-induced oscillations cannot safely be assumed to bear any simple relation- 
ship to flutter as it is known in the aircraft field, for they arise most often in a single 
degree of freedom and in an air stream which is not only unsteady but often mani- 
festly periodic in nature. An attempt has been made‘'’ to demonstrate that such 
an instability could occur with a cascade in incompressible potential flow, but only 
with such assumptions as to render the conclusions inapplicable in practice. 


Generally there is, even in the presence of significant aerodynamic forcing, no 
serious oscillation until the blading is driven beyond a critical combination of 
incidence and Mach number. In the possible operating range of a typical axial 
compressor there may be two or more such regions, one corresponding to stalling 
and another to choking flutter. Since the periodicity of the field in which the cascade 
is situated will have both incidence and Mach number components it may often 
produce a triggering effect, but one which need have no bearing on the subsequent 
oscillatory stress. 
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When a multi-stage compressor is throttled the incidence on its blading 
increases, but at a different rate in each stage. Thus, long before there is any one- 
dimensional instability of the general stream through the compressor some individual 
stage will stall, usually not as a whole, but in patches or cells which may move 
around the annulus at some fraction of the shaft speed’. The effect is of philo- 
sophical importance because it constitutes a marked distinction between the 
behaviour of infinite and finite but unbounded cascades. In practice it means that 
stalling flutter will often occur in regular but discrete bursts, while the blading in 
cascades other than that which is stalling will be subjected to what might with 
licence be described as periodically recurring, sharp-edged gusts. 


However, although each blade in the affected cascade may flutter only in bursts, 
it will be oscillating continuously as a result of the mechanical forcing transmitted 
through the disc or casing from other blades in the cascade. In the face of such 
physical complexity there is little hope of studying in detail the variation with time 
of the gas load on any one blade, so that the practical solution must lie in the 
adoption of a panacea either for the elimination of the operational circumstances in 
which flutter is possible, or for the mitigation of its effects when they occur. One 
possibility is the incorporation in the blade or its root fixing of artificially enhanced 
mechanical damping, and the principal object of this paper is a demonstration of 
the fact that the effects of such damping can be accurately assessed without reference 
to the aerodynamic mechanism of the flutter. 


Before proceeding it is desirable to consider in a little more detail the differ- 
ences between an aerodynamically forced oscillation and one which is self-induced. 
One practical distinction is that a damped blade may have quite a “ broadly tuned ” 
resonance curve and so could be forced (mechanically or aerodynamically) to 
oscillate over a wide frequency range, while in flutter the same blade would 
exhibit only its natural frequency. A more important distinction concerns the 
harmonic or overtone frequencies, for a blade can be forced not only in the region 
of its fundamental frequency but also in the region of one or more of its harmonics. 
For example, in the presence of the forcing wave form of Fig. 4 and Table II a 
typical blade would exhibit simultaneously its fundamental frequency, certainly its 
first, and perhaps also its second overtone—just like a musical instrument. In 
flutter, however, the energy will be extracted from the air stream principally in the 
region of the blade tip and the co-existence of an overtone mode would clearly 
upset the whole process. Thus, even when the periodic gas bending moment varies 
temporally in accordance with Fig. 4 and Table II, the blade will perform only in 
its fundamental mode and only at its fundamental frequency. This point is vital 
to the whole analysis and is why the results do not apply in instances of mechanically 
or aerodynamically forced vibration. 


Finally, to those skilled in the art of flutter as it is understood in the aircraft 
field, the mathematical approach and in particular the initial equations of this paper 
may seem a little strange. The answer might be said to lie in the fact that, whatever 
may be the form of its aerodynamic excitation, the motion of a fluttering cascade 
blade is known from observation to be very nearly sinusoidal. 
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NOTATION 


A, _ lifting moment Fourier coefficient 
B, lifting moment Fourier coefficient 
c “velocity” damping coefficient 
F Coulomb damping coefficient 
k _ effective stiffness 
m effective mass 
n term number in series 
p height of forcing pulse 
P, _ lifting moment Fourier coefficient 
t time 
x displacement of blade tip 
X  semi-amplitude of blade tip 
X positive part-amplitude of blade tip 
X_ negative part-amplitude of blade tip 
y “velocity” damping factor 
I’ magnification factor 
8 logarithmic decrement 
4 “F” vector phase lead 
vy term number in series 
p__ relative natural frequency 
 aphase angle 
a phase angle 
natural frequency in rad. /sec. 


2. The Basic Theory 
2.1. THE Exact THEORY 


Consider a single cascade blade which is oscillating in its fundamental flexural 
mode in an air stream of unknown state with a tip amplitude of (X,— X_) about its 
rest position as in Fig. 1 so that, whatever the aerodynamic circumstances, the 
periodic gas bending moment which it experiences is defined at any time ¢ by 
the series : — 


oo 
Periodic lifting moment= A, cosnwt+ B,sinnot. (1) 
2,3 3 


n=1, n=1,2, 


The period of this motion is 2*/w and the Fourier coefficients are in general 
functions of X;, X_ and the velecity and state of the air stream. 
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(MAGNITUDE OF DEFLECTION EXAGGERATED) 


Fig. 1. 
The physical system. 


In general, it will be necessary to assume for the mechanical damping of the 
blade and its root, and so on, both “ velocity ” and Coulomb components, the latter 
being represented for the moment by the square-topped wave: — 


Coulomb damping moment = >> sinvet ©& 
v=1,3,5 


which has a semi-amplitude of F. It should be observed that this assumption is 
valid only if, for dx/dt=0 at t=0, we also have dx/dt=0 at t=/w, and to this 
extent the equation of motion must anticipate its own solution. Should it transpire 
that there are regions in which the solution will not satisfy these conditions simultane- 
ously, the Fourier series representing the Coulomb damping moment will become 
generally indeterminate. However, for any numerical example, an exact solution 
could still be obtained by an iterative process involving the unknown coefficients. 


With these assumptions the equation of motion of the blade tip is then 


tin > Sinvet 
3 =1, 


d’x dx 
a y=1,3,5 


mie 


Ma 


to 


A, cos(nwt+9)+ 
,3 a= 


2 1 


where m, c and k are respectively the effective mass, “ velocity” damping coefficient 
and stiffness of the blade. The constant phase angle, ¢, is introduced to satisfy for 
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simplicity the boundary condition dx/dt=0 when t=O, and the coefficients m, c, k 
and F must be cast in moment form to satisfy both equation (3) and the physical 
system of Fig. 1. As shown in Appendix I, the “steady state” solution of 
equation (3) is then 


1 P, 
4F 1 1 —v2p? 
(1 — v2 p?)? +4 sin f— 
2p 
Tk y =1,3, s (1 — v2p?)? +4 COS (4) 
where P,= +(A,?+B,7)! 
tan (¢— Yn) = 252) A, (5) 
and 9 is given by 
S An—2nypB, .. x (I~ nip?) By 
4F 1 -—v? 6 


y=1,3,5 (1 —v*p?)? 


However, since this value of ¢ satisfies the condition dx/dt=0 when t=0, we can 
write for the positive part of the amplitude of motion 


1 P,, COS Wn 4F 


In this solution the non-dimensional parameters y and p are defined by 


and p= (8) 


2 (mk) (k/my 
That is, y is the ratio of the actual “velocity” damping coefficient to that which 
would effect critical damping in the absence of any Coulomb force, and p is the 
ratio of the frequency of oscillation to the natural frequency of the undamped 
blade. Unless they are actually negligible, the aerodynamic and centrifugal con- 
tributions to k and the aerodynamic contribution to m must naturally be considered 
in the estimation of the natural frequency. 


2.2. THE PHYSICAL SIGNIFICANCE OF THE EXACT THEORY 


Since the ultimate aim is the use in practice of an equivalent harmonic 
excitation, the most important feature of the exact solution is the rate of convergence 
of its series, for only if this is rapid will any reasonable approximation be achieved. 
There is thus a temptation to assume p=1, but a fairly rapid convergence is still 
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Fig. 2. 
The exact solution. 


obtained when p is given a value corresponding to oscillation at the natural 
frequency of the fully damped blade. This assumption, which is probably the 
closest approach that can be made to physical reality, is effected by writing 


so that the solution becomes a function only of the parameters 


A, B F 

In this form it can be applied directly and with equal facility to either flexural or 
torsional motion, for, apart from yy which is non-dimensional, the remaining 
parameters are effectively those steady deflections or twists which would result at 
the blade tip from the application of the moments A,, B, and F. 


Apart from the fact that its “F” vector anticipates a subsequent approximation 
(see Section 2.3), an interpretation of the solution for a typical instance is displayed 
in Fig. 2. However, if the offending vector were replaced by a series of rapidly 
diminishing ones, each with an appropriate phase angle and speed of rotation, the 
exact solution would be represented. From Fig. 2 the physical significance of » and 
vn is clear because, unless the condition dx/dt=0 when t=0 can be satisfied, the 
square wave series for the Coulomb damping will be invalid and must be replaced 
by a generally indeterminate function of wt. However, in numerical work ¢ is 
essentially the phase angle between the datum of the Fourier series representing 
the aerodynamic excitation and that representing the motion of the blade tip. That 
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it is therefore not in general the physical angle by which the peak of the exciting 
force leads the time of maximum tip deflection is more clearly shown in the example 
of Fig. 4. 


The example of Fig. 2 is one in which the excitation comprises three significant 
harmonics, but in order to study the influence of the Coulomb damping it is desir- 
able firstly to consider a case of simple harmonic forcing. In such circumstances y, 
would shrink to about two or three degrees and, A remaining unaltered, the point 
marked M would lie precisely on the motion datum line and we should have 


approximately 
k cos wf, 
where I" is a single-valued function of y. The physical effect of the damping is 


therefore to reduce the amplitude of the motion without distorting its simple 
harmonic wave form which is invariant with respect to both y and F. 


However, in the example displayed it is evident that, as the value of F is 
increased for a given value of y, the vector OM will in general attain zero magnitude 
at one or more values of wt while, at t=0, X, is still real and finite. The physical 
meaning of this is that, as shown by den Hartog? for simple harmonic forcing off 
resonance, the motion will proceed, but with one or more stops during each cycle. 
This, of course, destroys the validity of the series assumed for the Coulomb damping 
and so sets a limit to the relative values of y and F which can be considered by 
the simple theory. 


Finally, it appears from the diagram of Fig. 2 that, with the values of 9 and y, 
arranged as shown to give dx/dt=0 at t=O, it will not in general be possible to 
satisfy the condition dx/dt=0 at t=x/ and, moreover, that X, and X_ will not be 
equal. However, provided that there is only “velocity” damping the solution is 
still exact, but in the presence of Coulomb damping the equation of motion itself 
is wrong and the general theory becomes invalid. Nevertheless, as is shown in 
Section 3.3, a simple approximate method happens in practice to yield fairly 
accurate results. 


2.3. A SIMPLIFYING APPROXIMATION 


Consider now the magnitude of the Coulomb damping term of equation (4), 
namely, 
4F 3 1 —v?p? 
Tk y=1.3,5 (1—v?p?)? y*p? 
4F 2p 
Tk v=1,3,5 —v*p*)?+ 4v? 


Coulomb term= sin vot — 


5 COS 


The rapidity with which these series converge depends upon the value of y and is 
greatest for y—>0, so that we should consider first an example in which y is as 
great as it is likely to be in practice, say 0:2. 
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The magnification factor. 


We then have numerically, 
Coulomb term = — {2-538 cos (wt + 5°50’) +0-044 cos (3 t— 80°45’) + 
+0009 cos } 
o—2:538 cos (w t+ 5°50’). 


In making this approximation the error incurred at t=0 is only 0-2 per cent. 
and, for any value of t, the maximum error in the length of the “F” vector as shown 
in Fig. 2 is 1-6 per cent. Since this latter inaccuracy cannot affect the amplitude 
directly—but only through the validity of the Coulomb damping series of 
equation (2) when the amplitude is necessarily negligible—and the 0-2 per cent. error 
at t=0 will affect ¢ and y,, but little, the exact solution can safely be replaced for 
all practical purposes by 


P 


4F 
/ 
(1 — cos(nwt+v,)—-I cos(wt+A) (10) 


nk 


where ¢ is given by 


(1—n?p?) A, —2nypB, sin 2nypA,+ (1—n’p?) 


(1 —n*p?)’ + 4n? n=1,2,8 (1—n?p?)?+4n? 


cos¢=T sin A 
n=1,2, 


(11) 
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TABLE I 
[T’ AND JA OF THE APPROXIMATE SOLUTION 
y 0-01 0-025 0:05 0-10 0-15 0-20 
50:00 20-00 10-01 5-019 3-362 2:°538 
A 0 0°43’ 1°26’ 2°53’ 4°20’ 5°50’ 


and the functions I‘ and A have the values given in Table I and plotted in Fig. 3. 
As is more evident in Section 4, |" is merely the so-called magnification factor of 
orthodox forced vibration theory. 


3. A Limiting Example 
3.1. PRINCIPAL FEATURES OF THE EXAMPLE 


From the arguments of Sections 2.2 and 2.3 it follows that, when excited at 
resonance by a simple harmonic force, a blade will itself execute very nearly simple 
harmonic motion for all values of F at which the motion is significant and for all 
values of y up to about 0:2. In practice, any perturbation from simple harmonic 
motion must therefore be due to the higher harmonic content of the excitation and 
the form of the solution indicates that all such perturbations will increase in 
magnitude with y. The limits of approximation to the exact theory must therefore 
be found by studying an example in which the excitation is wildly non-harmonic and 
the damping coefficients are varied throughout the practical range. The chosen 
example is illustrated in Fig. 4 and, when it is expressed as a Fourier series, signi- 
ficant coefficients of the exciting moment are as listed in Table II. 


+ 
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TABLE II 
MOMENT COEFFICIENTS FOR THE LIMITING EXAMPLE 
n 1 2 3 4 5 6 7 8 100 
A,/p| —0:225 -—0159  -—0-075 0 +0:045 +0053 +0032 0 0 
B,/p| -0159 —07159 -0:109 -—0053 -0013 0 —0-006 
P,,/p 0°243 0°225 0-196 0-159 0-118 0-075 O 0:006 


Although physically impossible, this example can be thought of as an idealised 
case of stalling flutter in which the instantaneous lift, which is invariant with respect 
to x except at the points wt=0 and «/4, is known exactly and the problem is to 
estimate the amplitude and wave form of the resultant blade motion. Near the 
positive peak displacement there occurs, on the upper surface of the blade, a 
sudden separation of the boundary layer which lasts for a time «/(4) and produces 
a reduction of magnitude p in the root gas bending moment. Oscillation would 
occur at maximum amplitude when, due to the inertia of the separating air stream, 
the lift pulse lags behind the positive peak displacement by an amount ¢ which is 
precisely such as to satisfy the theory of Section 2. It will be observed that the 
first coefficient of the series accounts for only a small fraction of the total exciting 
moment and, whereas this may generally be so in stalling flutter, the example has 
been chosen with the intention that it shall in this respect be as bad as any case 
likely to arise in practice. Although relatively large, the hundredth harmonic is 
nevertheless unimportant, since what appears to matter is the distribution of 
magnitude amongst the first three or four harmonics only. 


3.2. SOLUTION FOR y=0-01 
Since for pure “ velocity” damping we have, by definition, 


the value of 0-01 for y yields a logarithmic decrement of 0-0628 as compared with 
the value for a typical compressor blade without root damping of about 0-001. The 
exact solution for F=0 (see Appendix II) is then as in Table III, so that we have 
for the positive part of the amplitude, 


12-17 — 0-00220 — 0:00139 —0:00080 — 0:00040 — 0-:00020 —0-00004...... ) 


12-177 (1 ~0-0050). 


TABLE Ill 
Exact SOLUTION OF EXAMPLE FOR y=0:01 AND F=0 
n 1 2 3 4 5 6 7 
12°17 0-0749 0-0245 00106 0:0049 0-0021 0:0007 
{(1—n?p?)? + 4n?y?p?}4 
Vn +0°50’ —110°59’ —133°46’ —156°24’ —178°24’ +158°29’ +135°36’ 
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Fig. 5. 
The motion without Coulomb damping for y=0-2. 


Moreover, although dx/dt~0 when t=7/, the error involved in overlooking this 
fact is clearly negligible because of the small value of ¥,, so with considerable 
accuracy we can write 


12177 (1-0-0014) 


and deduce for the amplitude of motion 


(X, —X_)=24-34 (1—0:0018)=24:307 


Despite its relatively small magnitude, the fundamental harmonic of the 
excitation thus gives rise to over 99-8 per cent. of the amplitude of motion, and the 
distortion from simple harmonic motion at the blade tip is likely to be negligible 


for all practical purposes. From equation (10) we can thus write with excellent 
approximation, 


(X,—X_)rao 
z k cos wf, 


whence the Coulomb damping moment required to reduce the amplitude to zero 
with the excitation unaltered is defined by F=0-191 p. 


3.3. SOLUTION FOR y=0-20 


For F=0 this amount of “ velocity” damping gives a logarithmic decrement 
of 1-282 and may reasonably be considered as beyond anything which will be 
necessary or perhaps even possible in future practice. The solution is defined by 
Table IV and, because of the relatively slow convergence of the series (see 
Appendix II), the wave form is distorted as shown in Fig. 5. 
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TABLE IV 
Exact SOLUTION OF EXAMPLE FOR y =0:20 AND F=0 
n 1 2 3 4 5 6 ij 
P 
0618 0:0763 0:0252 00110 0-00511 0-00223 0-00075 
{( p?)? +4n?+?p?} 4 


Vn +25°14’ —77°42’ —106°50’ —131°S4’ —155°42’ —179°7’ +157°23’ 


The boundary conditions for which equation (2) is valid are no longer satisfied, 
dx/dt being zero for ot=0 and + 155°, so that in the presence of Coulomb damping 
an exact solution can be found only by means of a rather lengthy computation. 
However, for F=0 the values of Table IV are exact and the amplitude of motion 
is therefore given by 


xX, =0:5597 (1 +.0-0291 — 0-0131 — 0-0131 — 0-0083 —0-0040—0-0012 +...... ) 
=0:559F (1 ~0-0107) 
X_=-0-618 (1 +0-0455) 


and (X,~X_)=1-177 (1 +0-0188)= 1-555 2. 
Thus, even under such disadvantageous conditions, the fundamental harmonic of 
the excitation would still account for over 98 per cent. of the periodic gas bending 
stress at the blade root. 


The error involved in the equation of motion when Coulomb damping is added 
is illustrated in Fig. 6, where it is assumed that the damping increment AF is too 
small seriously to affect either 9, ¥, or the amplitude of motion. If equation (10) 
be assumed valid, the calculated dissipation of energy due to any amount of 
Coulomb damping clearly cannot exceed that obtaining in fact and, for the small 
amount AF, we can write 


Actual dissipation per cycle = 2AF (X,—X_) 
Calculated dissipation per cycle= AF {[X, — X-]+ 


calculated dissipation 
actual dissipation 


so that =0-910. 

Now, if we make the gross assumption that they are valid, equations (10), (11) 
and (5) give the value of F required to stop the motion with the excitation unaltered 
as 0:158p and the resultant change in » and y¥, as about +6°9°. The consequent 
change in the sum of the P, vectors of Fig. 2 would thus be small, so that the 
difference between the Coulomb damping actually required to stop the motion and 
its calculated value of 0-158p would be due principally to the appearance in the 
series replacing that of equation (2) of an appreciable second harmonic term. 
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Fig. 6. 
The Coulomb damping inaccuracy. 


However, this term would produce such a small vector in the field of Fig. 2 that the 
exact critical value of F could reasonably be expected to lie some ten or at the most 
twenty per cent. below 0-158 p. 


3.4. THE PRACTICAL EFFECTS 

Suppose that in practice a blade with a value for y of 0-01 and no Coulomb 
damping is still found to exhibit a vibration with a tip semi-amplitude of, say, 
0-010 inches. If it is desired to reduce or eliminate this amplitude by an increase in 
damping, the form of the invariant excitation being unknown, there is no reasonable 
alternative to the assumption that the latter is harmonic. This being assumed, 
equation (10) reduces to 


Whence P/k is found to be 0-0020 inches. 


Now the results of calculations based on equation (13) will vary in accuracy 
as the unknown excitation varies in harmonic content and, since the example of 
Section 3.1 is as bad in this respect as anything likely to be encountered in practice, 
the maximum possible errors of such calculations can be found directly by com- 
parison with the results of Sections 3.2 and 3.3. Such a process yields the 
results of Table V. 


It thus appears that, in all practical instances of self-induced oscillation in one 
degree of freedom, estimations of the variation of amplitude with damping can 
safely be performed on the assumption that the unknown excitation is of simple 
harmonic form. The errors resulting from this superficially gross assumption are 
not only acceptable but in most instances actually negligible and, whenever the 
object is a reduction in amplitude, the error appears to be “on the safe side.” 
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TABLE V 
COMPARISON OF APPROXIMATE AND EXACT THEORY 


Circumstances | Result from Exact Theory | Result from Equation (13) Error 
{ Amplitude = 0-020 in. Amplitude = 0-020 in. Datum 
{ | i 157 micro-in. “i 157 micro-in. less than 1% 
{ Amplitude =950 micro-in. | Amplitude=1,016 micro-in.| +7% 
= F F 
=120to 130 micro-in. | — =157 micro-in. , 
Amplitude=0 | k o 130 micro-in k 157 micro-in +20 to 30% 


4. The Approximate Formula 


Should an aerofoil or blade be observed to oscillate aeroelastically in a single 
degree of freedom with a finite amplitude and a frequency defined by w rad./sec., it 
can be assumed to be excited by a periodic lifting moment of the form 


Excitation =P cos (wf + ¢) (14) 


where the appropriate phase angle, ¢, can for most practical purposes be taken as 
90° and P is in general a function of the blade tip amplitude and of the velocity and 
state of the air stream. Provided only that the variation of P with the semi- 
amplitude, X, is known, the variation of X with damping can then be estimated 
from the equation 


. 


1 


1 


The entities P/k and F/k are respectively those static deflections of the blade 
tip which would result from the steady application of P and the Coulomb frictional 
moment F, k being effectively the stiffness of the blade in its fundamental mode. 
The non-dimensional parameter y is the ratio of the actual “velocity” damping 
coefficient to that which would effect critical damping with F=0, 8 is the logarithmic 
decrement corresponding to y, and values of the magnification factor I" are given 
in Table I and Fig. 3. 


Provided that y does not exceed about 0:2, this value giving a logarithmic 
decrement of 1-282, the errors involved in the use of equation (15) will in practice 
be negligible, typical worst values being given in Table V. Moreover, by a simple 
change of co-ordinates, equation (15) can be applied with equal facility and accuracy 
to both flexural and torsional oscillations. 
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5. Concluding Remarks 


The formal conclusion of this investigation is stated in Section 4 and, whereas 
it is valid for self-induced aeroelastic oscillations, it must not be applied to instances 
of simple forced vibration, even when the forcing is purely aerodynamic. This 
is because of the basic postulate that the aerofoil or blade deflects in its fundamental 
mode and that this deflection is changed only in magnitude and not in form by the 
presence in the excitation of predominant higher harmonics. 


Before the method can be applied to the estimation of amplitudes, a know- 
ledge of the variation of P with X is of course essential, but this in no way affects 
the validity and accuracy of equation (15). Indeed, whereas throughout the 
argument it is implied that the method primarily concerns stalling or choking 
conditions where the amplitude of motion is finite, there is no reason why equation 
(15) should not be used to define the critical conditions for the onset in one degree 
of freedom of the diverging oscillation of a damped cascade blade, should such a 
phenomenon be discovered in practice (see, for example, Ref. 1). In such an 
instance, however, not only the centrifugal but also the aerodynamic contribution 
to stiffness might have to be taken into account. 
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Appendix I 
SOLUTION OF THE EQUATION OF MOTION 


As given in Section 2.1, the equation of motion of the blade tip is 


co ao 
dt? dt n=1,2,3 n=1,2,3 ® y=1,3,59 
(3) 
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Now the general solution of this equation comprises two parts, of which the first is an 
exponentially decaying term and the second is that which is known in forced vibration 
theory as the “steady forced vibration” solution. Only the latter is of interest here 
and, being a particular integral of equation (3), it is obtained by writing 


co oo co 
J, cosnwt+ K,sinnwt+SL,sinvot+ & M,cosvut (16) 


n=1,2,3 n=1,2,3 v=1,3,5 v=1,3,5 


On differentiating equation (16), substitution in equation (3) and a comparison of 
coefficients gives 


(k — mn?w?)J, +ncwK,=A, 
(k — mn’w?)K, ncwJ,=B, 


ign 
(k — my?w?)L, = ) 


(k — my?w?)M, +yvcwL, =0 


whence 
-1n?p?)A, — 2nypB, 
+ (1- n®p*)? + 4n?y2p? 


x 


K A, +(k—n’mw’*)B, 2nypA,+(1 - n?p?) B, 1 


(k — n?mw?)? + = n*p?)? + 4n?¥?p? k 
L k-v?mw? 1 vp? 1 
4F cw 4F 2yp 1 


However, it is convenient in practice to combine the first two terms of the solution 
vectorially so that equation (16) becomes, 


n=1,2, 


K 
(J,? +-K,,?)! cos( nt tan") +3 Lisinvwt+ M,cosyvwt . : (18) 
3 


n v=1,3,5 v=1,3,5 


and, on making this addition, the final “ steady state” solution of the equation of motion 
is found to be 


i Pe 


{(1 n°p?)? + 4n*y?p?}! Cos (nwt+Y,_)+ 
where P,, = +(A,?+B,?) 
tan“? 
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OSCILLATION OF COMPRESSOR BLADES 


an Only ¢@ now remains to be evaluated and, by definition, this is done by applying to the 
on solution the boundary condition dx/dt=0 when t=0, whence 
re 
A, 2nypA,+(l — n?p?) B, 
as (1 — n2p?)? + 4n?y2p? (1 n2p?)? + 4n?2p? COs 
6) 4F ©& 1 - v?p? 6 
(i - v2p?)? + 4y2y7p? (6) 
of 
Appendix II 
NUMERICAL SOLUTION FOR THE LIMITING EXAMPLE 
I) The system under consideration is illustrated in Fig. 4 and the Fourier coefficients 


of its exciting moment are given in Table II of Section 3.1. For the purposes of 
computation it is convenient to write 


1-—n*p?=G and 2nyp=H, 


and to know that p=(1-y*)}. 


EXAMPLE 1. y=001. 


TABLE VI 
(1) n 1 6 7 
(2) G? +H? 0:0004 9-002 64 22°59 576 1225 2304 
(3) GA, 0 0°477 0-600 0 — 1-080 £855 — 1-536 
(4) HB, —0-00186  —0-00636 -—0:01086 -—0-01272 -—0-01090 -0-00636 -—0-00182 
(5) (3) —(4) 0:00186 0-483 0-611 001272 —1:069 — 1-849 — 1-534 
(6) HA, —0:00450 -—0-00636 -—0-00450 0 0:00450 0:00636 0:00448 
(7) GB, 0 0-477 1-448 2°385 2°616 1-855 0°624 
(8) (6)+(7) —0:00450 1-444 2°385 2°621 1-861 0:629 
tan (¢—¥,,)=(8)/(5) | —2°419 09751 2°363 187-5 — 2-452 -1:0065 -0-4100 
(g—¥,) — 67°32’ 44°17’ 67°4 89°42’ 134° 49’ 157°42’ 
GA,—HB 
= ©) n 4-650 0:1073 0:0286 0:0002 -—0:0093 —0-0091 —0°0047 
G?+H? (2) 
H 
= ®) n| —11:25 0°1046 0:0677 0:0424 0:0228 0:0091 0-0019 
G? +H? (2) 


From Table VI, since A, = B, =0, we have 


- HB _HA,+GB, _ 


G?+ H? n=1,2,3 G?+H? 


whence, for F =0, equation (6) gives 
4-736 sin ¢ + 11:00 cos ¢=0, 
so that the numerical value of @ is — 66°42’. 


As is evident in Fig. 4, this superficially peculiar value arises as a-result of the 
position in time of the excitation series datum and has no physical importance. 
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In reality, the physical angle by which the peak of the exciting pulse leads the maximum 
deflection of the blade is given by 


Physical angle of lead = x — (66°42’) — (semi-pulse width) 


66°42 
= 90°48’, 


The vectorial solution, as illustrated in Fig. 2, would then be defined by the 
coefficients of Table II of Section 3.1 and the following phase angles :— 


n 1 2 3 4 5 6 7 
Wn 0°50’ —110°59’ — 133°46’ —156°24 — 178° 53’ 158°29’ 135°36’ 


EXAMPLE 2. y=0°20. 


The corresponding results for this amount of damping are as follows:— ” 
TABLE VII 
(1) a 1 2 3 4 5 6 7 
(2) G?+H? 0°15526 8-6795 59°751 208-67 532°8 1131-8 2126:2 
(3) GA, — 0:00900 0°45156 0:57300 0 —1:03500 —1:7787 —1-4733 
(4) HB, —0°03646 -—0°12458 -—0-21268 —0°24931 —0-21364 —0-12466 —0-03566 
(5) (3)—(4) 0:02746 0°57614 0°78568 0:24931 -—0°82136 -—1:6540 
(6) HA, —0-08820 —0-08813 0 0:08820 0°12466 0:08778 
(7) GB, —0-00372 0°45156 1-3828 2°2832 2°5070 1:7787 0°59852 
(8) (6)+(7) —0°09192 0:32698 1:2947 2:2832 2°5952 1:9034 0°68630 
tan (¢—¥,)=(8)/(5) | —3-3474 0:56754 1-64787 9-1581 —3-1596  —1:15079 —0-47739 
(¢—Y¥,) — 73°22’ 29°34’ 58°45’ 83°46’ 107°34’ 130°59’ 154°29' 
Pt 0-1769 0°1328 0:0395 0:0048 0:0077 0:0088 —0-0047 
G+H? 
HA,+GB (8) 
n @) n 0:5920 0:0754 0650 0:0438 0:0244 1 


Thus, for F=0, 
0-3326 sin + 0-3712 cos ¢=0, 


whence ¢= — 48°8’ and the physical angle by which the exciting pulse leads the maxi- 
mum blade displacement is 109°22’. As before, the solution is then defined by the 
coefficients of Table II of Section 3.1 and the following phase angles : — 


n 1 2 3 4 5 6 of 
25°14 —77° 42’ — 106°50’ — 131°54’ — 155°42’ —179°7’ 
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The Drag of Non-Planar Thickness 
Distributions in Supersonic Flow 


E. W. GRAHAM, B. J. BEANE and R. M. LICHER 
(Douglas Aircraft Company, Inc., Santa Monica, California) 


SUMMARY: For many conventional aircraft, the thickness of the wing and 
fuselage can be represented approximately by source distributions in the 
plane of the wing and on the axis of the fuselage. Such aircraft may be 
considered as having essentially planar thickness distributions. However, 
missiles with cruciform wings, biplane arrangements, “ring” wings and so 
on, require non-planar source distributions to represent the wing and fuselage 
thickness. For this reason, the term “ spatial thickness distribution ” is used. 


To simplify the investigation of spatial thickness distributions, a 
singularity representing an element of volume is introduced. It is shown that 
the optimum distribution of such elements in a prescribed space gives rise to 
a minimum wave drag value consistent with that obtained for a Sears- 
Haack optimum body of revolution. 


It is also shown that the interference of lift and side force elements, 
with the optimum distribution of volume elements alone, makes it possible 
to reduce the drag of a given volume below the Sears-Haack “ minimum.” 
Such interference effects (among lift force, side force, and volume elements) 
apparently account for the nearly zero drag of finite aspect ratio Busemann 
biplanes, and the zero drag of Ferri’s axially symmetric central body plus 

outer shell. 


1. Introduction 


The distribution of lift on an aircraft in supersonic flight can be represented by 
the use of elementary horseshoe vortices (lifting elements). The thickness of the 
wing, fuselage and other components can be represented by source distributions. 
These singularities satisfy the linearised equation for supersonic flow. 


Many conventional aircraft can be represented fairly well by distributions of 
sources and lifting elements in a single horizontal plane. In such cases there is no 
interference between the sources and lifting elements, and the problems of reducing 
drag due to lift and drag due to thickness can be studied separately. 


However, the investigation of more general types of aircraft requires the use 
of singularity distributions which cannot be confined to a single horizontal plane. 
The term “spatial distribution” is introduced to cover this broader class. In 


Received September 1954 


May 1955 99 


n 

i- 

e 


E. W. GRAHAM, B. J. BEANE AND R. M. LICHER 


addition to lifting elements and sources it is necessary to consider side force 
elements, which are simply lifting elements rotated through 90° so that the forces 
are directed to one side instead of vertically. 


In the most general case, these singularities (sources, lifting elements and side 
force elements) have interference drag terms, so that the drag due to thickness, for 
example, cannot be studied independently of the lift and side forces on the aircraft. 


In this paper, the thickness drag is considered first and the interference effects 
caused by lift and side force are then discussed. 


NOTATION 


a,b — semi-minor and semi-major axes of an ellipse 
c length of a cylinder, or length of body, or wing chord 
Cy drag coefficient 
d___ distance between source and sink 
D drag force 
f() 
f, intensity of singularity distribution 
F() 
K aconstant 
length of body 
l(x) | chordwise radial load distribution on shell 
m — mass of fluid introduced per unit time 
M free stream Mach number 
q dynamic pressure 


co-ordinate normal to streamwise direction, or radial co-ordinate on 
circular plan form wing 


R _ semi-span of circular plan form wing 
S cross-sectional area 
t thickness 
U, _ free stream velocity 
volume 


co-ordinate in streamwise direction (X = ax/R) 


= 1-(r/RY 
= /(M?-1) 


# angular inclination of the trace of the Mach planes on the xy-plane, 
relative to the y-axis 


V 
x 
y co-ordinate in spanwise direction (Y = by/R) 
a 
B 
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€ co-ordinate in the streamwise direction or co-ordinate perpendicular 
to trace of Mach planes in the xy-plane 


» co-ordinate parallel to trace of Mach planes in the xy-plane 


Rasa) 


potential 


¢, potential of a unit source 
potential of a unit dipole 


6 polar co-ordinate on wing plan form or angular co-ordinate from 
the xy-plane 


Pp, free stream density 
Optimum 


()min minimum 


Dashes denote differentiation 


2. The Singularity Representing an Element of Volume 


The investigation of lift distributions is simplified by the use of a singularity 
which represents an element of lift. This singularity is the elementary horseshoe 
vortex. The intensity of lift corresponds to the strength of the singularity and the 
location of the lift force is identical with that of the bound vortex. The study of 
volume (or thickness) distributions is similarly simplified by identifying the singu- 
larity which corresponds to an element of volume. 


Consider a source and sink of equal strength and located on the same stream- 
wise line. In each unit of time a certain quantity of fluid is introduced into the flow 
pattern by the source and the same quantity is removed by the sink. The volume 
occupied by the fluid flowing from source to sink depends on the strength of the 
source and sink and the distance between them, and also depends on the velocity 
and density of the fluid flowing from source to sink. However, if the volume is to 
be considered a linear function of the strength of the singularities, then the mean 
value of density times velocity must be unaffected by the perturbation velocities 
created by the source and sink. This means that in a linearised treatment of the 
problem the fluid flowing from source to sink may be considered to have free stream 
density and velocity. 


Let m= mass of fluid introduced per unit time 
d= distance between source and sink 
Py = free stream density 


U,= free stream velocity. 


Then the volume occupied by the fluid is a ‘ 
0 0 
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Since the volume is proportional to md, doubling the intensity of source and 
sink and halving the distance between them should produce a shorter, but thicker, 
volume of the same magnitude. This suggests proceeding to the limiting case (as in 
incompressible flow) where the source and sink are combined in a dipole with axis 
in the free stream direction. This singularity should represent an element of volume, 
although the fineness ratio of the element is zero. 


The potential for a unit source at (¢, 0) in supersonic flow is 


-1 


where B= /(M?-1), 


x and € are co-ordinates in the streamwise direction and r is radial distance from 
the x-axis. 


Differentiating with respect to x gives 


+(x-& 
2x[(x— £)?— 


+ dy, 
where ¢, is the potential for the unit dipole or an element of volume equal to 1/U,. 


3. The Distribution of Volume Elements 


~~ - 


Fig. 1. 


For a distribution of volume elements along the €-axis with intensity f(é), 
starting at £=0, the potential is 


2m 


Integration by parts gives 


The first term in the expression for the potential is infinite, and apparently 
corresponds to the “ roughness” of the body, which is an assembly of blunt elements 
(see Fig. 1). 
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The smoothly faired body (indicated by broken lines) is all that we are 
concerned with, and this creates the finite part of the potential. This finite part is 
also the potential for a source distribution of intensity equal to +f’ (£). This source 
distribution can be used to construct a body of revolution extending from —//2 to 
+1/2. The volume of this body, computed from the surface shape obtained from 
slender-body theory, is equal to 


+1/2 


{ 


—1/2 


This is consistent with the total volume corresponding to the sum of the 
volume elements. 


Similarly, a planar distribution of volume elements should give a wing with a 
thickness distribution whose total volume computed from the surface shape is equal 
to the sum of the volume elements. 


The concept of the volume element is not particularly useful in studying bodies 
of revolution and planar wings, since these cases are relatively simple. However, 
the use of the volume element does help to clarify problems involving more general 
spatial distributions of thickness. 


The points to be emphasised are that fixing the sum of the volume elements 
fixes the total volume and fixing the distribution of volume elements determines the 
drag. It is therefore possible to study the drag of a distribution of volume elements 
without calculating the exact shape of the corresponding body. This is analogous 
to the fact that the drag of a distribution of lifting elements can be studied without 
calculating the twist and camber of the corresponding wing surfaces. 


4. The Drag of Volume Distributions on a Streamwise Line 


A body of revolution may be constructed from a distribution of volume 
elements along a streamwise line, or from the equivalent distribution of sources. 
The body constructed from volume elements is an “ infinitely rough” body and has 
infinite drag. However, discarding the infinite part of the potential leaves a 
“smooth” body (with finite drag) which is equivalent in every respect to the body 
created by a source distribution. 


If f (x) is the intensity of the volume element distribution for a body of revolu- 
tion of length /, then the drag is given by"? 
+1/2 41/2 


D= | og |x, x, | dx, dx, 


-—1/2 


To maintain constant total volume according to linearised theory, 
+1/2 
[ f (x) dx=constant. 
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The body shape giving minimum drag for a given length and volume has been 
determined by Sears’ and Haack’”’ independently. The corresponding f(x) (which 
is proportional to the cross-sectional area) is given by 


+1/2 
| 
fon | | | — |. 
Thus, the optimum distribution of volume elements along the axis corresponds 


to the cube of an elliptical distribution. (For lifting elements the optimum distri- 
bution is elliptical). 


The value of the minimum drag is 


+1/2 2 


f (x) dx 
Duin= 7 (4) U, (l/2)° 


= 2G) Cal 


5. The Drag of Volume Distributions in Space 


Fig. 2. 


If the volume elements are not confined to a single streamwise line, then the 
drag contributions at different angles, 6, on Hayes’s cylindrical control surface? 
are not necessarily the same. For any one angle, 0, the drag is given by 


41/2 41/2 
39 = | f’ (x,,9 f” (x2, 9) log | x, — x2 | dx, 


—1/2 


Here f(x, #) is determined by the use of “ Mach planes” for the angle 6. All 
the volume elements intercepted by any one “ Mach plane” are transferred (in the 
plane) to the streamwise axis. The resulting distribution along the axis is f(x, 4). 
The problem of finding the minimum drag contribution at the one angle @ is then 
equivalent to the Sears-Haack problem. If f(x, 9) corresponds to the cube of an 
elliptical distribution for a// 6, then the total drag is a minimum, and the drag 
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contribution at each 6 is a minimum and corresponds to that of an equivalent Sears- 
Haack body. 


It is not always possible to get such a perfect minimum. For example, if the 
space within which thickness may be distributed corresponds to that which can 
be utilised in a “ planar” wing of rectangular plan form, then it is not possible for 
the drag contributions at all angles 6 to be made a minimum simultaneously. How- 
ever, if the available space corresponds to that occupied by a “planar” wing of 
elliptical plan form, it is possible to get a “perfect” minimum. The optimum 
distribution of thickness for such a wing was first determined by R. T. Jones". 
(In Appendix I of this paper, the optimum thickness distribution for a wing of 
elliptical plan form and constant volume is determined by a slightly different 
procedure). Jones found that for minimum drag with constant volume the aerofoil 
sections should have constant curvature in the streamwise direction. The resulting 
thickness distribution, viewed from any angle on Hayes’s cylindrical control surface, 
gives an equivalent body of revolution having the Sears-Haack optimum shape. 
However, these various equivalent bodies, corresponding to different angles on the 
control surface, have different lengths, so there is no one Sears-Haack body which 
corresponds exactly to the elliptical wing. 


Fig. 3. 


If we consider the optimum distribution of thickness within a space which has 
axial symmetry about a streamwise axis, then it is possible for all the equivalent 
Sears-Haack bodies to have the same length. For example, consider that a double 
Mach cone bounds the space within which thickness is to be distributed (Fig. 2). 
The Sears-Haack body placed on the axis has the same drag contribution at every 
angle on the cylindrical control surface and, of course, the “equivalent” body of 
revolution for any angle @ is identical with the real body. However, a “ring” wing 
(which carries no radial forces) plus a central body of revolution (Fig 3) can be 
designed to have exactly the same drag as the Sears-Haack body (see Appendix 
II(A)). The equivalent bodies of revolution are all identical with the Sears- 
Haack body. 


In obtaining the minimum drag values just discussed, only volume elements 
were used. However, it will be shown in the next section that lift and side force 
elements have interference drag with the optimum distribution of volume elements. 
Consequently, the drag for a given volume can be reduced below the Sears- 
Haack value. 
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6. Interference Between Lift and Side Force Elements and 
an Optimum Distribution of Volume Elements 


Consider a Sears-Haack body placed on the axis of a double Mach cone, and 
place a lifting element as shown in Fig. 4. The interference drag between body and 
lifting element is composed of two parts, the effect of the body nose on the lifting 
element and the effect of the lifting element on the tail of the body. 


The nose of the body corresponds to a source distribution and produces an 
upwash velocity at the lifting element. This causes negative drag. The lifting 
element produces a positive pressure at the tail of the body. This also causes 
negative drag so the total interference drag is negative. (This argument, of course, 
applies not only to the Sears-Haack shape but to other shapes also). 


LIFT ELEMENT TRAILING VORTICES 


The total drag of the combination is equal to the drag of the Sears-Haack body 
alone plus the drag of the lifting element alone plus the interference drag. The 
drag of the lifting element alone is proportional to the square of the lift it carries. 
However, the interference drag is proportional to the first power of the lift on the 
element and to the first power of the strength of those sources and sinks in 
the body which are affected by interference. The lift carried by the element can, 
therefore, always be made small enough so that the drag of the element alone is less 
(in absolute magnitude) than the interference drag. Thus, the total drag of the 
combination can be made less than the drag of the Sears-Haack body alone. 


This suggests placing elements of lift and side force in a ring about the Sears- 
Haack body, and so arranged that the force on each element is directed radially 
outward from the body. This process may be used to construct a central body plus 
cylindrical shell which has zero drag (Appendix II(B)). Such a system has been 
investigated by Ferri‘; its two-dimensional analogue is the Busemann biplane. 


It therefore appears that the optimum distribution in space of volume elements 
alone yields minimum drag values consistent with the Sears-Haack values. How- 
ever, the optimum distribution of volume elements plus lifting and side force 
elements should give zero drag for any total volume. 
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Appendix I 
OpTIMUM THICKNESS DISTRIBUTION FOR ELLIPTICAL PLAN FORM 


It is desired to find the optimum thickness distribution for a planar wing of elliptic 
plan form and given volume; this problem was first solved by R. T. Jones’. A 
geometrically simpler problem, which will be examined first, is to find the optimum 
thickness distribution for a circular wing of given volume. The method of Hayes", 
in which the drag is evaluated by summing increments of drag at each angular station 
around a cylindrical control surface far away from the body, will be used. For the 
total drag to be a minimum, the increment of drag at each angular station should also 
be a minimum. Hayes’s method substitutes for the actual planar distribution, as viewed 
at each angular station, an “ equivalent’ thickness distribution concentrated along an 
axis parallel to the free stream direction. It is known that the optimum thickness 
distribution for such an equivalent body is that of the Sears-Haack body":?). For 
each angular station the equivalent body is constructed by finding the cross-sectional 
area of the wing intercepted by a series of planes inclined at the Mach angle to the 
free stream direction. The area cut out by a given Mach plane is then concentrated, 
at the intersection of the Mach plane and longitudinal (streamwise) axis, such that the 
volume of the equivalent body is equal to the volume of the actual wing. 


If the thickness distribution of the circular plan form is axially symmetric, then 
the equivalent bodies at each angular station will have the same shape (although 
different “ fineness ratios”) due to symmetry. If t(r) is the thickness distribution which 
is to be the optimum for a given volume V, then 


R 
0 


where R is the radius of the circular wing and r, @ are polar co-ordinates from the wing 
centre (Fig. 5). The area cut out at each point along the £-axis by planes normal to 
that axis is 


+ (R? R 


The equivalent lineal distribution along the x-axis is 
R 
t(r)rdr 


S(x)=2 cos (3) 
COS 
with { S(x)dx=V. 
—Rsecp 
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Up 
x 
RSEC 
Fig. 5. 
For minimum drag, this distribution should be“ 
3 
x 2 
Thus the integral equation to be solved for t(r).is 
R 
t(r)rdr 
Rsec u (r? — x? cos? 1)’ @) 


COB 


where K is a constant dependent upon the given wing volume. By a suitable trans- 
formation of co-ordinates, equation (4) may be written in the form 


K t (a) da 


x 2 
where g=l1- ( ) 
Rsec p- 


r 2 
R 
Equation (5) is called Abel’s equation and its solution is well known, cf. Ref. 6. 


The solution to equation (5) is 


| 


and substitution of this in equation (1) determines K; then 


(a) 
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y 


+b 
4 
U, 
Fig. 6. 


Equation (6) thus gives the distribution of thickness which will result in minimum drag 
for the circular plan form wing of given volume. 


To apply the circular plan form solution to the original problem of finding the 
optimum thickness for an unyawed elliptic plan form, we make the foliowing change 
of co-ordinates : — 

X=ax/R 
(7) 


Y=by/R 


The circular wing is then transformed into an elliptic wing whose equation is 


It can be verified that the thickness distribution 


obtained from equation (6) through the transformation of equation (7) is the optimum 
for this more general case; that is, the equivalent lineal distribution for equation (8) 
with a set of Mach planes inclined at the angle « as shown in Fig. 6 is 


where l= J (a? + b? tan? 1). 


Since equation (9) represents a Sears-Haack body, the thickness given by equation (8) 
is optimum for the unyawed elliptic wing.* 


Determination of the total drag in this optimum case involves an integration of the 
drag increments from these Sears-Haack bodies as seen at each angular reference station. 
If the reference station is at an angle 6 from the horizontal, then the Mach planes cut 
the elliptic plan form at an angle » defined as'*? 


*Within the accuracy of linearised theory, this thickness distribution corresponds to constant 
curvature in the streamwise direction. 
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and the total drag is D= | — do. 
The increment of drag at each reference station is‘) 


log|x — | dx dé 


4n 


x 8 
dD_ -q | 
do dx? dé 
0 


4qV? 
and the total drag for the optimum thickness distribution of equation (8) is 
4qV? | (M?-1)+2a?/b? 
Defining t=t, (1 X?/a? Y?/b?) 
and D=Cpgqrzab, 
M? — 1) + 2a?/b?]t,? 


[(M?— 1) a 


This result agrees with that given by Jones“). 


Appendix II 


(A) RING-WING AND CENTRAL BODY OF REVOLUTION COMBINATION HAVING DRAG OF 
SEARS-HAACK BoDy 


Consider a ring-wing plus a central body of revolution contained within the volume 
bounded by a double Mach cone. Because of the axial symmetry of this particular 
system, the equivalent body of revolution is independent of the angle 6 on the cylindrical 
control surface. In this case, if the local radial force on the wing is everywhere zero, 
the drag of the equivalent body of revolution is, according to Hayes’s formula, identical 
to the drag of the original system. Thus, a ring-wing (which carries no radial force) 
plus a central body of revolution will have exactly the same drag as a Sears-Haack 
body if the equivalent body of revolution is a Sears-Haack body. 


To design such a system, we may select any smooth, slender profile for the ring- 
wing and compute the cross-sectional areas cut from this wing by a set of parallel Mach 
planes. These areas must then be subtracted from the cross-sectional areas which would 
be cut from a central Sears-Haack body by the corresponding Mach planes. The 
resulting area difference defines the area distribution (in the Mach planes) of the correct 
central body. (This area must be projected normal to the flow direction to obtain 
the cross-sectional area of the central body defined in the usual way). This body, 
together with the ring-wing originally selected, is an optimum distribution of thickness 
within the double Mach cone volume. 
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DRAG OF NON-PLANAR THICKNESS DISTRIBUTIONS 


He 
RING WING (t/C= 0-06) 
= CAMBER REQUIRED BUT NOT 
SHOWN (RADIAL FORCE=0) 
- 


ENTRAL BODY 


OF REVOLUTION 


Up 
EQUIVALENT BODY OF 
REVOLUTION (SEARS-HAACK, 
FINENESS RATIO=6) 


t 1 


1 


Fig. 7. 
Cross-sectional view of ring-wing and central body. (An optimum distribution of thickness 
within the double Mach cone volume). 


o/c 


As an example, consider a ring-wing with thickness distribution corresponding to a 
bi-parabolic arc profile. The camber necessary for zero local radial force need not be 
determined, since it does not affect the shape of the central body. Assume that the wing 
is six per cent. thick and located half-way between the axis and the apex of the volume. 
If the central body of revolution is designed so that the equivalent Sears-Haack body 
is of fineness ratio 6, the resulting shape of the central body of revolution is as shown 
in Fig. 7. 


(B) RING-WING AND CENTRAL BoDy OF REVOLUTION HAVING ZERO DRAG 


The theoretical minimum drag value for a distribution of thickness elements that 
has no interference with lift or side force elements is the drag of a Sears-Haack body. 
Part (A) of this Appendix illustrates a system which attains this minimum drag value 
for thickness elements contained within a double Mach cone volume. However, it was 
pointed out previously that interferences between thickness distributions and distributions 
of lift or side force may provide negative drag contributions which reduce the theoretical 
minimum wave drag of a system to zero. For the double Mach cone volume, the 
remaining portion of this Appendix illustrates a central body of revolution which, 
together with a certain distribution of radial forces on a cylindrical shell, has zero wave 
drag. The method used here to design such a system makes use of certain equivalences 
between sources and line distributions of elementary vortex shells.* (These equivalences 
will be discussed in a later paper). 


*A vortex shell is defined to be-a ciosed vortex line in a plane normal to the flow direction. 
The present application assumes infinitesimal rectangular vortex shells distributed uniformly 
along a streamwise line. 
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CYLINDER 


FILLED 
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NG | 
/ 4 


DOUBLE MACH 
CONE VOLUME 


Fig. 8. 


Consider an axially symmetric, continuous distribution of sources filling a cylindrical 
space contained within the double Mach cone volume (Fig. 8). A set of parallel Mach 
planes cutting this source distribution determines an equivalent lineal source distribution 
according to the method of Hayes’. For convenience, this equivalent lineal source 
distribution will be denoted by F(x). Because of the axial symmetry of this case, 
F(x) is independent of the angle 6 on the distant control surface. A central body of 
revolution which is represented by the negative of F(x) will just cancel the velocities 
induced at the distant control surface by the original sources. The drag of the combined 
system is then zero. The remaining step is to relate the original source distribution to 
a distribution of radial forces around the boundary of the cylindrical space. It can be 
shown that a source and a sink of equal strength, lying on the same line parallel to the 
flow direction, have exactly the same effect at the distant control surface as a line of 
constant-strength elementary vortex shells connecting the source and the sink. If such 
vortex shells are considered to replace the original source distribution, the forces on 
adjoining shells inside the cylinder are seen to cancel one another, while the forces on 
the outer sides of the shells next to the boundary of the cylindrical space determine the 
radial force. A cylindrical shell having this radial load distribution plus a central body 
of revolution which corresponds to the source distribution — F(x) constitute a system 
having zero drag. 


As an example, suppose that a cylinder within a double Mach cone volume is 
considered to contain a source distribution which varies linearly with axial distance 
but is independent of radial distance. That is, the source strength per unit area inside 


the cylinder is 
f=-4,(1-2%), 
Cc 


where x is measured from the leading edge of the cylinder, c is the cylinder length, and 
f, is the strength of the sources at the rear face of the cylinder. The equivalent lineal 
source strength corresponding to this original distribution is given by 
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t U(x) 


RADIAL FORCE DISTRIBUTION 
ON CYLINDRICAL SHELL 


O-5}+ 


CROSS-SECTIONAL SHAPE 
DOUBLE MACH CONE OF CENTRAL BODY OF 
REVOLUTION WHICH CANCELS 
(AT DISTANT CONTROL SUR- 
FACE) DISTURBANCE DUE TO 
M2/2 RADIAL FORCE 


j 
-05 0-5 1-5 


Fig. 9. 
Central body of revolution and radial force distribution having zero wave drag. 


and -2(1-22)]; BR<x<c- BR, 


where R is the cylinder radius. Now, if the negative of this source distribution is 
assumed to represent a body of revolution, then within the accuracy of slender-body 
theory the area distribution of the central body is 


S(x)= - r [ F (x) dx. 
— BR 


For illustration, the dimensions of the cylinder are assumed to be such that 
BR/c=1/2; that is, the radius is half the distance between the axis and the apex of the 
double Mach cone volume. The shape of the central body of revolution which cancels 
the effect of the original source distribution for this case is shown in Fig. 9. The 
distribution of radial force which can replace the original linearly varying source 
distribution is 
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Supersonic Flow Past Slender Bodies With 
Discontinuous Profile Slope 


L. E. FRAENKEL and H. PORTNOY 
(Imperial College) 


SumMMarRY: Ward’s slender-body theory is extended to derive first approxi- 
mations to the external forces on slender bodies of general cross section 
with discontinuous profile slope. Two classes of body are considered: 
bodies whose profile (typified by the local radius) is continuous between the 
nose and base, and certain bodies whose profile is discontinuous, such as 
bodies with annular or side air intakes and wing-bodies on which the wing 
has an unswept leading edge. (Where air intakes are concerned, it is 
assumed that they are sharp-edged and that there is no “spillage” of the 
internal flow). 


The following conclusions apply to the former class of bodies. The 
variation of drag with Mach number is found to depend only on the discon- 
tinuities in the longitudinal rate of change of the cross-sectional area, and is 
thus independent of cross-sectional shape. The drag itself is unchanged if 
the direction of the flow is reversed. The expressions for lift and moment 
assume the same forms as for smooth pointed bodies, the lift depending only 
on conditions at the base of the body. 


The general theory is applied to winged bodies of revolution with an 
unswept wing leading edge: the results bear a marked resemblance to those 
obtained by Ward. The results for wings alone are seen to be applicable, 

with one modification, to subsonic as well as to supersonic speeds. 


1. Introduction 

In this paper we extend Ward’s slender-body theory"? to derive first approxi- 
mations to the aerodynamic forces in steady supersonic flow on slender bodies of 
general cross section whose profiles have discontinuous slope. (By the word 
“ profile” we mean the boundary of a section of the body by a plane through its 
longitudinal axis). The discontinuities must lie approximately in planes normal to 
the direction of the undisturbed flow. Lighthill‘? first treated this problem for the 
particular case of bodies of revolution by introducing some results of the quasi- 
cylinder theory. In Ref. 3 bodies of elliptic cross section were considered, and it 
was shown that for a first approximation to the forces the introduction of a quasi- 
cylinder theory was not absolutely necessary; a modified form of Ward’s simple 
velocity potential was used in a particular manner to obtain the results. It is this 
approach which is generalised in the present paper, and it is found that the 
method may also be used for certain bodies whose actual profile is discontinuous. 
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The physical effect of the discontinuities is that in short regions behind them 
the flow is roughly analogous to that on a wing with a sharp supersonic leading 
edge, and the pressure is of a greater magnitude there than it is on smooth, pointed 
bodies. Consequently these regions, although short, make contributions to the 
aerodynamic forces which are of the order of the total forces. The mathematical 
effect on the slender-body approximation has been discussed in Ref. 4 and elsewhere: 
this solution now provides only an asymptotic approximation, valid when the 
distance downstream from a discontinuity is large compared with the thickness 
of the body, and it fails completely in the important regions of relatively large 
pressure behind the discontinuities. This difficulty is overcome by relating the 
integral over the body of the 0¢/0s term in the pressure to the values of the 
potential ¢ itself at points immediately upstream and far downstream of a discon- 


tinuity, where the approximation is valid. 


NOTATION 
a,(s) coefficients in the slender-body potential 
b(s) wing semi-span 
b,(s) function in the slender-body potential 
B J/(M?-1), where M is free stream Mach number 
C,C;, | body contour in a plane s=constant, s=e;+ 
D drag force 
e; values of s at which R, dv/ds are discontinuous (i=0, 1, 2, . . . m) 
F,G complex lateral force X¥ + iY as given by ¢, ¢ 
1 length of body 
1, value of s at maximum span of wing 
M.,M, moments about x- and y-axes 
M. complex moment M,+iM, 
(r,6,s) cylindrical co-ordinates 
R(6,s) value of r on body surface 
R,, Re _ partial differentials of R with respect to s and @ 
R,(s) radius of body of revolution 
real part of 
S(s) cross-sectional area of body 
t maximum thickness of body 
U _ free stream velocity 
w(z,s) complex function such that 2w=¢@ 
(x,y,s) Cartesian co-ordinates 
X,Y _ forces in positive x-, y-directions 
z,Z, complex variables x + iy, (x+ 2s)+iy 
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Z,(s) complex co-ordinate of centre of area of cross section 
a angle of incidence 
B(s) effective wing semi-span, (b*— R,”)/b 
V. _ two-dimensional gradient operator in plane s=constant 
vy normal outward from body in plane s=constant 
dv/ds streamwise body slope, RR, (R? + 
p _ density in the free stream 
circumferential variable, dr =(R?+ Ro’)! dé 
yg exact linearised perturbation potential 
@ asymptotic slender-body potential 
see equation (3) 


2. The Slender-body Solution 

Let (x, y,s) and (r,6,s) be Cartesian and cylindrical co-ordinates such that an 
undisturbed stream at infinity, of velocity U, flows in the direction of increasing s, 
and such that the origin of the frame is at the centre of area of the nose cross section 
of the body. Let the complete velocity potential of the field be ®=U (s+-¢). 


Consider a body subject to the restrictions specified by Ward", except that its 
streamwise slope dv/ds and its profile may be discontinuous at stations s=e; +O (¢’), 
(i=0, 1,...m), which must be spaced well apart (e:,,-—e;,=O(1)). We take e,=0, 
€m=l, and specify that / is O(1). The following conventions are also used : — 


fix. 
The linearised equation for the velocity potential is 


v2 

where V,, denotes the two-dimensional gradient operator in a plane s=constant; and 
the boundary condition on a slender body is 


We assume the existence of a perturbation potential ¢ which satisfies equations (1) 
and (2) exactly, and which also satisfies the appropriate conditions at infinity and 
on the wave fronts emanating from the discontinuities. 
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The asymptotic, slender-body solution #, valid near the body for s—e;» O (?), is 


where #, is harmonic in (r, 9), satisfies the boundary condition (2), and is such that 


(s)logr]—>0 as r—> oo; and 


Here a Stieltjes integral has replaced the Riemann integral of Ref. 1. 9%, may be 
expanded in series outside a circular cylinder enclosing the body : — 


[ (s) log z+ a,(s)z-" |, r>Rrax, 
n=1 (5) 
where a, (s)= S’ (s). 


The coefficients a, (s) are O(t"+*) and are in general discontinuous at s=e;. 


It may be shown"? that if we write -=9+¢,, then there exists some finite 
positive number K such that 


for S—e,<Kt, 9, is o (tos 
oe 


while for s—e,>Kt, 9, is o +). 
™ —~ Oi 


The arguments leading to this last result include an assumption which is supported 
by known particular cases. 


3. The Drag Force 


The drag integral is first expressed in terms of ¢, and cast into a form which 
permits substitution of the asymptotic potential ¢. The drag of a segment of the 
body between s=e; and s=e,,, may be written 


% 


The O-term results from the fact that in the initial region s—e; <= O(), where 
02/0s is O(t), the error of @¢/0s is O(t’), because ¢ is only a solution of the 
linearised equation. 
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We introduce the total derivative of the integral around the body contour of 
¢(d¢/dv). Let C denote the body contour in a plane s=constant, C’ the projection 
of C on the plane s+és=constant, and C” the body contour in the plane 
s+és=constant. Then* 


aim § [2 

= lim, +¢ 2) Os Ov +? Hes 
(7) 
where we have used Green’s theorem for the first integral, and the fact that C’ is 
geometrically the same as C, but is displaced a distance 4s, in the second. Proceed- 


ing to the limit 5s —> 0, and substituting for V,?¢ and 0¢/0v from equations (1) and 
(2) respectively gives 


Inserting this result into (6), we now have 


dv 

In this expression the contribution of 2 B?2 (0*¢/0s?)(dv/ds) is at most O(t*). In 
the second integral we substitute 


9, e,—)+O (#* (log? 


since ¢ (as opposed to ¢) is continuous, and in the remaining terms we may substitute 
¢ =: the error is still O(t°) (see the Appendix). 


The integral over the body surface may now be transformed, by the two- 
dimensional form of Green’s theorem, into an integral over the surface of a circular 
cylinder r=r, > Rmax, ¢:<5<é;,,, (Fig. 1), on which the series expansion of ¢ 
converges. We have 


26 9 ar r,d0—2 _.0)(7, 


Or ds 
0 
=2| Spas 1247 dS— os ~ dr. (10) 


* The derivative 0*:/@vds appearing in (7) may be defined as (v. V,) (0¢/ 0s). 
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Fig. 1. 


Inserting this expression into the double integral of (9), integrating with respect to s, 
and transforming the resulting integrals of (V,¢)? over A, and A, into integrals 
around the boundaries of these areas, we obtain 


+2 | a dS — r,d6+O(t*). 


(11) 


Substituting ¢=9,+ 5, in the first two integrals of this expression, and using the 
series expansion of ¢ in the others, we now find 


This result in no way requires the profile function R (6, s) to be continuous at s=e;; 
it only requires the existence of a slender-body type of flow field immediately 
upstream of C;,, on s=e;—, r=R(6,e;+). We shall, however, only sum the drag 
of the m segments for bodies for which R (6, s) is continuous in the range O0< s</. 
Performing this summation, and inserting the value of b, (s) given by equation (4), 
we obtain, after some reduction, 


D 


1 2 


g=0— 
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where ¥Y denotes the “finite part” of the double Stieltjes integral; that is, the 
points s=c=e, are excluded from the area of integration (see Ref. 2). Two 
general conclusions may be drawn from this result : — 


(i) The variation of drag with Mach number depends only on the discon- 
tinuities in the function S’ (s); it is independent of cross-sectional shape and 
of the area distribution between discontinuities. 


(ii) The-drag is unchanged if the direction of the stream is reversed: this is 
proved as follows. In reversed flow the signs of S’(s), ¢, and dv/ds 
are all changed, and values at s=e;+ are replaced by values at s=e,— 
Hence it is only necessary to show that 


That this is in fact an identity may be shown by applying Green’s theorem 
to the area between C; and the circle r=r,, s=e,, and using thé series 
expansion of ¢ on the circle. 


4. The Lateral Forces and Moments 


The complex lateral force due to a segment of the body is 


={as$ [- 2% _ |ide+ (14) 


Fi (X + i+: 
4pU* 4pU* 


Proceeding in a manner similar to that of equation (7) we find that 


=> +50; 5 st) dz |. 


where we have used (dv/ds) dr=(0¢/0v)dz =(02/0x) dy —(0¢/0y)dx. Hence (14) 
becomes 


= 


We may again substitute ¢ ~@ into this expression, taking the value at s=e,— in 
the first term. The inner part of the double integral is then 


ow\? _ 
7 (3) dz, 
120 


fag (52 +152) . (15) 


e; 
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which vanishes, since 0w/0z is a regular function of Z outside C, and since the 
integral around an infinite contour vanishes. Taking the most general case where 
R (0@,e,-) A R (0, e,+) we find the lateral force up to some station s at a substantial 
distance downstream from a discontinuity. The difference between the integrals 
around the body contours immediately upstream and downstream of a discon- 
tinuity in R is transformed into an area integral by Stokes’ theorem, and there 
results 


d 
=4na,+2£ (Sz)-2 Hs-e)| (16) 


where H is the Heaviside unit function, and the last step has been taken directly 
from Ref. 1. The total complex lateral force is obtained by evaluating this 
expression at s=/-. 


To calculate the moments, let G(s) denote the complex lateral force up to s 
as given by ¢, and let F(s) denote that given by ¢. Then the total moment is 


M.=M,+iM,=i sG’ (s) ds + O(t*) 


0 


G(s) being a continuous function. It is permissible to substitute G (e)) ~ F (e;—), 
G (s) = F (s), and hence 


M.=i[ IF(-)- [F(@ds| 


0 


Thus for bodies on which R(@,s) is continuous in the range 0<s</, the 
expressions for lateral force and moment are unaffected by the discontinuities, and 
the lateral force still depends only on conditions at the base of the body. 


The results of Refs. 2, 3 and 5 can, of course, be recovered from equations (12), 
(13), (16) and (17) as particular cases. 


5. Winged Bodies of Revolution 

In this section the results are applied, in the manner of Ref. 1, to calculate the 
forces on winged bodies of revolution at incidence. The novelty of the present 
treatment lies in the fact that the body geometry and the wing plan form may now 
include discontinuities (Fig. 2). The wing has zero thickness and total span 25(s), 
and is mounted symmetrically on the body, the incidence of wing and body being 
defined by z,= —2s; the plan form of the configuration must be slender; and the 
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Fig. 2. 


discontinuities must be separated by distances of O(1). If s=/, denotes the station 
of maximum wing span, the radius of the body R,(s) (now measured from z,), 
must be constant for s>/,. In general it is necessary also to impose the restriction 


b(s) b’(s)=>R, (S)R,' (8), << (18) 


which ensures that the side edge of the wing is a “ leading” rather than a “trailing ” 
one with respect to the local flow, so that the vortices do not drift outboard of the 
wing, upstream of the maximum span. However, when R,’(s)=0 for s>e,, or 
when the wing alone is considered, this restriction is probably not necessary. In 
these latter cases the body does not deform the streamlines in the wing plane, and 
if b’(s) is negative (so that ],=e,) the unknown flow over the foremost part of the 
wing only affects the vorticity in the most outboard part of the wake, the width of 
this part being O (t’). 


The asymptotic potential of the wing-body is” 


R { a, log z, + b, (s)+ al (z.+ 0<s<l,, 
(19 


b* (s)—R,? (8) 


BG) (20) 


where Z,=z+es, B(s)= 


8 (s) might be called the “effective wing semi-span,” (see next paragraph). 


We insert this potential into equations (12), (16) and (17), noting that there is no 
coupling between the thickness and incidence terms in the resulting contour integrals, 
due to the anti-symmetry of the incidence terms. Writing D, for the drag of the 
body at zero incidence, we obtain 


Pate 
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28 {1,s«,)- [eas] }. (23) 


0 


The results for a wing alone may be of some interest. Such results are of course 
a particular case of those just given (D,=S=0, e,=0, B=b), and can also be 
obtained, for b’(s)>0, from those for a body of elliptic cross section®, by con- 
sidering a body of suitable profile and length /,, and letting the minor axis of the 
ellipse tend to zero. 


The results for the lift and moment of a flat slender wing are valid for subsonic 
speeds also, since essentially they depend on only two facts: (i) that ¢, provides an 
asymptotic approximation to ¢, (b, (s)=0 here), and (ii) that in the plane of the wing 
¢=0 upstream of the leading edge. Both these remain true for subsonic speeds. 


The drag result given here cannot be applied to subsonic speeds because the 
suction force at the unswept part of the leading edge cannot be calculated by the 
methods used. The subsonic drag can be found, however, by a well-known method. 
If we calculate it by considering the flow of momentum, as predicted by ¢, through 
a large circular cylinder surrounding the wing, the contributions of the upstream 
face and of the curved surface of the cylinder vanish in the limit. There remains 
the contribution of the downstream face or “Trefftz plane.” In this plane ¢ is 
harmonic in (r,@) and the drag is 


where the contour C,, is taken around the wake. Hence 


which is the same result as for a pointed wing. 


Most of these results for a wing alone were originally suggested by the work 
of Jones’, and were confirmed, for the case of a slender rectangular wing in 
supersonic flow, by Stewartson’. 
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Appendix 


In this Appendix we justify the substitution ¢=9@ in the term 
32 
7 
e; c 


of equation (9). The proofs for other, similar substitutions proceed along the same 
lines. For convenience we put e,=0, e;,,=/,. 

We introduced the error ¢, of the asymptotic potential in Section 2; we now 
establish the order of 0°¢,/d0vds on the body. Differentiating along a line s=constant 
locally, we have 


ds\ dvas dv? ds’ 


But ¢ and @ satisfy the same boundary condition; hence the left hand side of this 
equation vanishes. Applying the orders of @, given in Section 2 to 079, /dv?, we now have 


=O tog for s>Kt. 

Ss 
In the region s<= Kt, we may say that 079,/dvds is O(t) at most, since 0*¢/dvds and 
0*9/dvds are both O (t). Hence 


ly Rt 


Jas > ? ds log TOM dr+ 


0 Cc 0 Cc 


1 


2 
+ Jas [o+0(S tos! 


Kt Cc 
1, 


Kt 
4, { j at } 
Cc 


0 0 Kt 


log 


and the O-term is found to be O (f°). 
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The Buckling and Bending of Orthotropic 
Sandwich Panels With All Edges 
Simply-Supported 


J. R. ROBINSON, B.Sc., A.F.R.Ae.S. 
(Handley Page Limited) 


Summary: A formula is derived for the buckling in compression of an 
orthotropic sandwich panel, with all edges simply-supported. Two special 
types of symmetrical sandwiches, that of a sandwich having zero core flexural 
stiffnesses and of a corrugated core sandwich, are considered, and curves 
giving the critical buckling end load per inch are presented. An approximate 
theory is also given for the bending of sandwich panels under the action of 
an axial compressive load and a uniform lateral pressure. Curves are 
presented giving the maximum bending moments in a corrugated core 
sandwich panel. 


1. Introduction 


In the design of sandwich skins for aircraft wings, one of the important features 
is the buckling and bending of panels. The difference between sandwich theory 
and the normal plate theory is the low transverse shear stiffnesses which can no 
longer be neglected and must be included in all sandwich calculations. While most 
of the reports on sandwich construction deal with isotropic cores, this paper deals 
with an orthotropic core. 


The assumption is made that the sandwich behaves essentially as a plate, so 
that it can be defined by the following physical constants : — 
(i) the two flexural stiffnesses, D, and D,, 
(ii) the twisting stiffness, D.,, 
(iii) the two transverse shear stiffnesses, Do, and Dg,, 
and (iv) the two Poisson’s ratios, », and »,, defined in terms of the curvatures. 


The buckling and bending of an orthotropic sandwich panel, with all edges 
simply-supported will be considered. The mathematical theory is given in the 
Appendix, while in the main text the steps of the mathematical theory are explained 
and only the important equations are quoted. 


The numbers beside the equations in the main text are those used in 
the Appendix. 
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Fig. 1. 
Forces and moments acting on differential element 5x dy. 


A 


NOTATION 
x, y, Z rectangular co-ordinates; z is measured normal to the plane of 


the panel (see Fig. 1) 

w deflection of the middle surface of the panel, measured in the 
z-direction (in.) 

q__ lateral loading (1b./in.*) 


Q. internal shear force per unit width, acting in the z-direction 
in a cross section originally parallel to the yz-plane (Ib./in.) 


Q, internal shear force per unit width, acting in the z-direction 
in a cross section originally parallel to the xz-plane (Ib./in.) 


M,. __ internal bending moment per unit width, acting in a cross 
section originally parallel to the yz-plane and about a line 
parallel to the y-axis (lb.in./in.) 

M, __ internal bending moment per unit width, acting in a cross 
section originally parallel to the xz-plane and about a line 
parallel to the x-axis (Ib.in./in.) 

M., internal twisting moment per unit width, acting in a cross 
section originally parallel to the yz-plane or xz-plane 
(Ib. in. /in.) 


N. _ internal tensile force per unit width, acting normal to a cross 
section originally parallel to the yz-plane (Ib. /in.) 


N, _ internal tensile force per unit width, acting normal to a cross 
section originally parallel to the xz-plane (Ib./in.) 


N., shear per unit length in the plane of the panel and parallel to 
the yz-plane and the xz-plane (Ib. /in.) 


D,, D, _ flexural stiffnesses of the panel for bending about the y-axis 
(= bending moment (lb. in.) 


curvature 
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Rectangular plate with horizontal forces applied to the boundaries. 
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twisting stiffness of the panel for twisting about the x-axis 


( twisting moment Linch) (Ib. in.) 
twist 


shear stiffnesses of the panel (1b. in.) 


Poisson’s ratios for the panel (defined in terms of curvature) 


Poisson’s ratio for the face plates 


length and width respectively, of the rectangular panel (in.) 
(see Fig. 2) 


number of half waves in the x-direction 


number of half waves in the y-direction 


constants defined by equations (17)-(21) and (42) 


4n°*D, 
b? (1 — pe py) 


= P,/(4Dau) 


D,/D, 
Dgz/Dey 


D.,/D, 


a/b 
critical buckling end load per inch of the sandwich panel of 
infinite length (parallel to Ox), i.e. critical value of —N, 
maximum bending moment parameters defined by 
(Mz)max = 5°, (My)max= %yq 5? 
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2. Displacement and Equilibrium Equations for the Panel 


By using the small deflection theory, including the shear terms, a set of six 
differential equations can be obtained, three of which express the equilibrium of an 
infinitesimal plate element and the other three relate the curvatures and twist of the 
element to the forces and moments acting on the element. A detailed derivation 
of these equations is given in Ref. 1 and only a brief one is given in the Appendix. 
These six equations can be reduced to the following three equations by eliminating 
the bending and twisting moments. 


wD, \ D. 1 #0, 
1 Diy ft, D, \ 1Dy 
* U2 Day * * 2 Dye 
Pa D, 0? 1 Duy 
2 Doz oxdy 2 De, 0x’ Q, (9) 
ay +N, ay? +2Noy dxdy) . (10) 


The problem now is to obtain a solution of these differential equations which 
satisfies the appropriate boundary conditions for the problem being considered. 


3. Boundary Conditions 
The boundary conditions for a simply-supported panel are: — 


(i) no deflection normal to the plate, i.e. parallel to the z-axis, 
(ii) no bending moments about the lines of support, 
(iii) no displacement of the face plates along the lines of support. 


4. Rectangular Panel Subjected to an End Load per inch, N,, 
Parallel to the x-Axis and a Varying Lateral Loading 
The first system of external loads to be considered will be an end load per inch, 
N,, parallel to the x-axis and a lateral loading g of the form 


MTX 
4 = nn Sin sin =, 


where m and n are integers. 


128 The Aeronautical Quarterly 


8) 


9) 


0) 


ch 


rly 


BUCKLING AND BENDING OF ORTHOTROPIC SANDWICH PANELS 


The solution of the differential equations (8)-(10) for this system of external 
forces is of the form 


mee sin 
b 


W=Ann Sin 


0.= Bun cos” 


(12) 


. mn 
QO, =Ca. sin 


These equations satisfy the boundary conditions for a simply-supported panel, 
and the constants Amn, Bnn and Cy». can be obtained by substituting the equations 
(12) into equations (8)-(10). When these constants are known the deflections, 
internal forces, and moments of the plate can be obtained. 


5. Buckling of an Orthotropic Sandwich Panel for Compression 
Parallel to Ox 
The buckling load of the panel can now be obtained from the condition which 
makes the deflection w infinite. This gives the following expression for N. in terms 
of the physical constants of the sandwich, the aspect ratio 8 of the panel, and the 
integers m and n. 


k 


In this equation, m and n are the number of half waves that the sandwich 
buckles into along and across the panel. The critical form of buckling across the 
panel will be when n=1, but that along the panel will depend on the aspect ratio 8. 
By plotting —N, against 8 with various integers for m, the minimum values of 
—N, and hence the critical buckling end load per inch can be obtained. 


(32) 


6. Buckling in Compression of a Symmetrical Corrugated Core 
Sandwich Panel 
One type of sandwich in present use is the corrugated core sandwich. Algebraic 
expressions for the elastic constants for this type of sandwich are given in Ref. 3. 


It can easily be proved that the shear stiffness along the corrugations is many times 
that across the corrugations, so that sufficient accuracy is obtained by putting 
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Fig. 3. 
Compressive buckling end load per inch of a symmetrical corrugated core sandwich panel 
(all edges simply-supported). 
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Fig. 4. 
Compressive buckling end load per inch of an infinitely long corrugated core sandwich panel 
(all edges simply-supported). 
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Kez 
_ 
105 0-2 0-3 04 . 0-6 0-9 1-0 
Fig. 5. 


Critical value of 8/m for an infinitely long corrugated core sandwich pancl (all edges 
simply-supported). 


= equal to infinity. This is justified by Fig. 9. With the substitution, equation (32) 


then reduces to 


kat kp? +1) 


k 
_ BD 


P, 
ltr {1+ B +e) 

Figure 3 gives the buckling parameter —N,/P, plotted against the panel aspect 
ratio 8 for three values of r, when the longitudinal flexural stiffness is 1-25 times 
the transverse flexural stiffness, which is representative of a typical corrugated core 
sandwich. The critical buckling end load per inch, P.,, of an infinitely long panel 
of width b is obtained by minimising (37) with respect to 8/m and is plotted in 
Fig. 4. The corresponding critical half wavelengths divided by the panel width are 
plotted in Fig. 5. 


7. Buckling in Compression of a Symmetrical Sandwich Panel, 
Having Zero Core Flexural Stiffnesses 


The other type in general use is the honeycomb core sandwich. This comes in 
the class of sandwiches having zero core flexural stiffnesses, that is the flexural 
stiffnesses D, and D, are equal. For this set of sandwiches equation (32) reduces to 


m BY 
N.. Bom (38) 


P, 
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Fig. 6. 


Compressive buckling end load per inch of a symmetrical sandwich panel having zero core 
flexural stiffnesses (all edges simply-supported). 


Figure 6 gives the buckling parameter — N,/P, plotted against the panel aspect 
ratio 8 for various values of r, which depends on the ratio of transverse flexural 
stiffness to transverse shear stiffness for two values of ¢, the ratio of longitudinal to 
transverse shear stiffnesses. The critical buckling end load per inch, P.., of an 
infinitely long panel of width b is obtained by minimising (38) with respect to 
8/m and is plotted in Fig. 7. The critical value of 8/m is plotted against r in Fig. 8. 
It is possible in minimising the expression for —N, in equation (38) to obtain 
negative values for (8/m)*. In these cases, the minimum value of —N, is when 
8/m equals zero and depends only on the shear stiffnesses of the panel, and is 
plotted as a dotted line in Fig. 7. The values of P.. corresponding to very small 
values of 8/m are only of academic interest. This corresponds to the case of 
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Fig. 7. 


Compressive buckling end load per inch of an infintely long symmetrical sandwich panel having 
zero core flexural stiffnesses (all edges simply-supported). 
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Fig. 8. 


Critical value of 8/m for an infinitely long symmetrical sandwich panel having zero core 
flexural stiffnesses (all edges simply-supported). 
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Fig. 9. 
Variation of P.,./P, with < for an infinitely long symmetrical sandwich panel having zero core 
flexural stiffnesses (all edges simply-supported). 


anti-symmetrical wrinkling and the assumption that the sandwich can be regarded 
as a plate is no longer valid. For these conditions some of the other physical 
constants of the sandwich, particularly the modulus of elasticity of the core normal 
to the face plates and the flexural stiffnesses of the face plates alone, become 
important, while some of the physical constants used are no longer important. The 
wrinkling of sandwich panels has already been covered by several authors. 


Figure 9 gives the critical buckling end load per inch, P.,, of an infinitely long 
sandwich panel for a fixed value of r, which depends on the ratio of transverse 
flexural stiffness to transverse shear stiffness. From this it is seen that the maximum 
possible increase in P.. by raising the ratio of longitudinal to transverse shear stiff- 
nesses, ¢, from 10 to a higher value, is less than five per cent. This fact allowed the 
simplification, by taking ¢ equal to infinity, of the calculation of the buckling loads 
for a symmetrical corrugated core sandwich panel with corrugations running parallel 
to the applied load. 


For an isotropic core 1 
It has been shown in Ref. 5, that the formulae derived by several authors for the 


critical buckling end load per inch of an infinitely long sandwich panel with simply- 
supported edges and isotropic core reduce to this expression. 


8. Rectangular Panel Subjected to End Load Parallel to the x-axis 
and Uniform Lateral Loading q 


The uniform lateral loading can be represented by the double trigonometrical 
series 


MEX. 
q= = SiN 
m=1,3,5, b 
16 
mn 
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Fig. 10. 
Values of (M,),-»/. for k= 1-25, r=0-333, B=3-0, and 
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Fig. 11. 


Maximum bending moment, (M,),,,,, in a rectangular corrugated core sandwich panel 
subjected to end load per inch parallel to the x-axis and a uniform lateral loading q 
(panel aspect ratio=3-0), 
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Fig. 12. 
Maximum bending moment, (M,),,,,, in a rectangular corrugated core sandwich 
panel subjected to end load per inch parallel to the x-axis and a uniform lateral loading q 
(panel aspect ratio=3-0). 


Expressing the uniform loading in this manner allows the equations for the deflected 
surface, shear forces and bending moments to be obtained by taking the summation 
of the corresponding values previously obtained. In the Appendix expressions are 
derived for the bending moments M, and M, (equations (49) and (50) for a 
corrugated core panel and equations (54) and (55) for a panel having zero core 
stiffnesses). On the numerical side, only typical examples of the use of this method 
for corrugated core sandwich panels are given because of the large amount of 
computing work necessary in deriving each curve. Fig. 10 gives the bending 
moment about an axis parallel to the y-axis, (M.),-»,., at points along the centre 
line of a typical corrugated core panel of aspect ratio 8=3-0, for the case when r, 
which depends on the ratio of transverse flexural stiffness to transverse shear stiff- 
ness, equals 1/3 and k, the ratio of longitudinal to transverse flexural stiffness equals 
1-25, for various values of the end load parameter —N,/P,. For this case it is seen 
that the peak bending moment occurs at approximately one-sixth of the length of 
the panel from either end. The maximum values of the bending moments M, and 


M, can be expressed as 


(Mz) max = b? (52) 


(M,)max = ay qb? (53) 
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Fig. 13. 


Maximum bending moment, (M,),,,;, in a rectangular corrugated core sandwich panel subjected 
to end load per inch parallel to the x-axis and a uniform lateral loading q. 
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Fig. 14. 


Maximum bending moment, (M,),,,,, in a rectangular corrugated core 
sandwich panel subjected to end load per inch parallel to the x-axis 
and a uniform lateral loading q (r=0°4, —N,/P, =0-7). 
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where z, and z, are functions of —N,/P,, r, k and 8. The values of 2, and 2, are 
plotted in Figs. 11 and 12 for a typical corrugated core panel, where 8=3-0 and 
k=1-25. The effect of changing the aspect ratio 8 of the panel on z, only has been 
computed because of the work involved, and is shown in Figs. 13 and 14. 
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Appendix 
Mathematical Theory 


1. DISPLACEMENT AND EQUILIBRIUM EQUATIONS FOR THE PANEL 


Consider a small element of the panel of length 6x and width dy, cut from the 
unloaded panel by planes parallel to the xz- and the yz-planes. The internal forces 
which act on this small element are shown in Fig. 1. A detailed derivation of the 
displacement equations is given in Ref. 1 and only a brief one is given here. 


The curvature 0?w/dx? of the panel in the zx-plane is made up of three 
components : — 

(i) that due to the bending moment M.,. and equal to — M,,./D.,. 

(ii) that due to the bending moment M, and equal to + »,M,/D, (i.e. anticlastic 

curvature). 

(iii) that due to the variation of the shear force Q, and equal to (0Q,,/0x)/Do,. 

The differential equation for the curvature 0?w/0x? is obtained by summing these 
three components, and is 

M, M, 1 0Q, 


Similarly the differential equation for the curvature 0?w/dy? is 

M, M, 1 0Q, (2) 

2 

dy? D, De, 


The twist, 0*?w/dxdy, of the panel is also made up of three components : — 


(i) that due to the twisting moment, M,.,, and equal to M,,,/D,,. 
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(ii) that due to the variation of the shear force, Q,, and equal to 


11 20, 
2 Der dy 


(iii) that due to the variation of the shear force, Q,, and equal to 


112, 
2 Doy Ox 
The differential equation of the curvature 0?w/dxdy is obtained by summing these 
three components and is 


aw M, 1 100, 1 1 @Q, 


dxdy Dz, * 2 Dos 2 Do, Ox 


(3) 


The next step is to obtain the equilibrium equations. It will be assumed that the 
middle-plane stresses are constant throughout the panel. The equation for equilibrium 
of vertical forces is given at the top of page 305 of Ref. 2 as 


02M 
dx? Axdy Oy? 


2 2 2 


(4) 


The equation for equilibrium of moments about the y- and x-axes is obtained from 
equations (188) and (189) of Ref. 2 as 


_ OM, 

oy (5) 
_ OM, 


These six equations (1)-(6) will be reduced to 3 equations by eliminating M,, M, 
and M,,. 


Eliminating M, between (1) and (2) gives 


aw 1, 1 2Q, 


Substituting (3) and (7) in (5) gives 


Do; dy? + 2Do, dxdy 0x NOx? 
1 0?Q 1 
— fly Da ox y Oxdy 
D, D, 1 D, 0°Q 


1 My 1 _ 
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Similarly, from (1), (2), (3) and (6), 


Me Dy \ 1 Dy _ 
{33 Doz 1- Py Doz Oxdy 2 De Ox? 2, 9) 
and from (4), (5) and (6), 
ax oy Ox? oy? +2Noy 


By assuming general values for the middle-surface deflection w and the shear forces 
Q, and Q, which satisfy equations (8), (9), (10) and the relevant boundary conditions, 
it is possible to solve a number of bending and buckling problems of sandwich panels. 


2. BOUNDARY CONDITIONS 
The boundary conditions for a simply-supported edge parallel to the y-axis are 
(i) no deflection parallel to the z-axis, i.e. w=0. 
(ii) no bending moment about a y-axis, i.e. M,=0. 


(iii) either there is no y-displacement of points in the boundary, so that the shear 
angle Q,/Do,=0; or the only support is applied to the middle surface of the 
boundary when M,, =0. 


The edge supports for sandwich panels are usually such that there is no displacement 
of the face plates along the line of support. Hence the boundary conditions for the 
rectangular panel shown in Fig. 2 are 


(i) for the sides x=0 and x=a, 


w=0 
Q,/Dgy=9 
(ii) for the sides y=0 and y=b 
w=0 
Q,/Dg,=90 
3. RECTANGULAR PANEL SUBJECTED TO AN END LOAD PER INCH N,, PARALLEL 
mux 
TO THE x-AXIS AND A LATERAL LOADING OF THE FORM SiN 


The assumption will be made that the solution for the middle-surface deflection 
w and the shear angles Q,/D 9, and Q,/ Dg, exists in the form 


nt. 
w=A,,, Sin —— Si | 
a b 


Q,=C,, sin— — cos 
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These satisfy the boundary conditions (11) and (12). 
Substituting (13) in (8), (9) and (10) gives 


+ )+™ D, 


B mx? D, 1 1 n?x? D,, }- 
(1 My } 
nine Dy, } 


mnt? (1 Dry Be 


-1;=0 1 
{ b? 1 | 2 a? ( 5) 
ir bus mn mn =0 (16) 
le 
at Solving equations (14), (15) and (16) for Aina, Bina and C,,, gives 
1e 
Amn — Bus Cun — 1 
where 
Ww (= Dz _ min 
) m'n?x* D,D, 
a*b? (1 pe py) 
1) 1 Dz 1 D, m 1 n? 
1 D,D., mn? DD, (uzD,+ 
— + + Dow @ 1-pep, pop, (19) 
23 D 


4 2 2 D 2 


mr nt 
Zan = ma) Xn 21) 
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Hence the equations for the deflected surface and the shear forces are 


Z.. a b 
MTX nny 
The bending moment M,, is, from (7), 
D, [ 0 fow Q, ) 0 (ow Q )| 
—[ — - - — 2 
De, dy \0y Dey 
Substituting (22) in (23) gives 
D, [ mx 1 nt ] 
M, 1 Mer mn a Xin a By h2 mn bDoy x 


D m X. PPh Y 
1 a’ + My b? a Doz b Dey 


Amn... 
Similarly 
m? nY 

—“ sin sin pa (25) 

a b 


mn 


The expressions for W,,,, Xin» Ym, and Z,,, are simplified by the following 
substitutions : — 


4x*D 
r=P,/4Do, 
¢=Do,/ Dey 
n= D,,/D, 
B=a/b J 


Also from the generalisation of Maxwell’s reciprocal law 
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Substituting (26) and (27) in (18), (19) and (20) gives 


r 


k BF 

[3 ky Bt B (k - 1 + | + 

1. xm [ m (k-p,”)+ 

+ Pu kB +n? ( (29) 


2 2 
lp m n (k Pe \+ku, 


(30) 


Hence the deflected shape, shear forces, and bending moments can be obtained by 
substituting equations (28)-(30) in equations (22), (24), and (25). 


4. BUCKLING OF AN ORTHOTROPIC SANDWICH PANEL FOR COMPRESSION PARALLEL TO Ox 


The end load per inch, — N,, which will cause the panel to buckle can be obtained 
by putting the denominator of the expression for w (equation (22) ) equal to zero. 


ive. (== Wan ) 
a b a 


Substituting (28), (29) and (30) in (31) gives 


+n? (14 +1] 


r(k—-p,? 1 n> m? 


1 =2m?n? n (k—u,2)+k 


The minimum value of — N, given in (32) is now required. The critical form of buckling 
in the y-direction will be when n=1. 
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This reduces (32) to 


m* 
tha +25 k +1 (33) 


By plotting - N, against 8 with various integers for m, the critical buckling value of 
- N, of the panel is obtained. 


5. BUCKLING IN COMPRESSION OF A SYMMETRICAL CORRUGATED CORE SANDWICH PANEL 


For a symmetrical corrugated core sandwich panel with the corrugations Tunning 
parallel to Ox, the following substitution can be made in (33):— 


Dz,=D,/(1 +p) 
For other details of the Elastic Constants for a corrugated core sandwich see Ref. 3. 


From (26) and (34), n= Ita (35) 
Substituting (34) and (35) in (33) gives 
k-p?+ky+kp? 
(36) 
m? (k 
(= +(1 * = 


Normally for a corrugated core sandwich Dg, is many times greater than Do, and 
sufficient accuracy is obtained by putting < equal to infinity. This is justified by Fig. 9. 
With this substitution equation (36) reduces to 


+& 2m? k- +kut+kp? +1) 
4N,  kd+p) 
P, = (37) 
B 
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The buckling parameter - N,/P, is plotted against 8 in Fig. 3 for the case where 
k=1-25 and »=0-25. 


The buckling end load per inch P,,, of an infinitely long panel of width b is obtained 
by minimising (37) with respect to B/m. P., is plotted in Fig. 4, taking 1.=0-25. The 
corresponding critical half wavelengths divided by the panel width (i.e. 8/m) are plotted 
in Fig. 5. 


6. BUCKLING IN COMPRESSION OF A SYMMETRICAL SANDWICH PANEL HAVING ZERO CORE 
FLEXURAL STIFFNESSES 


For a symmetrical sandwich panel having equal flexural stiffnesses, k=1-0, and 
the substitutions defined by equations (34) still hold. These reduce equation (33) to 


1 f/m 3 
N 
2 
e B 26 B 


is (+) 


The buckling parameter, — N,/P,, is plotted against 8 for two values of ¢ in Fig. 6. 
P., for an infinitely long panel of width b is obtained by minimising (38) with respect 
to B/m and is plotted in Fig. 7, taking »=0-25. The critical value of B/m is plotted 
against r in Fig. 8. It is possible in minimising the expression for — N, in equation (38) 
to obtain negative values for (@/m)*. In these cases, B/m is taken as zero to obtain 
the critical values of P,., which will depend only on the shear stiffnesses of the panel. 


To show the effect of the ratio of the shear stiffnesses, «, on P,,, a curve has been 
drawn in Fig. 9 for the case r=0-5. From this curve it is seen that in increasing « from 
10 to infinity there is less than five per cent. change in the critical buckling end load per 
inch. This was the justification for taking « equal to infinity for the corrugated core 
sandwich panel. 


The buckling end load per inch for an isotropic core is obtained by putting «=1-0 
in equation (38) which gives 
1 A) 
tm 


q l+r(1+ 


For an infinitely long panel P,, is obtained by minimising (39) with respect to B/m 
and is given by 


(39) 


P, 


This formula is similar to that derived by several authors for the critical buckling 
end load per inch of a simply-supported sandwich panel, with an isotropic core (for 
example, see Ref. 4). 
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7. RECTANGULAR PANEL SUBJECTED TO END LOAD PARALLEL TO THE X-AXIS AND A 
UNIFORM LATERAL LOADING qg 


The uniform lateral loading can be represented by the double trigonometrical series 


MTX nk 
q= > Ann Sin —— sin ; (41) 
m=1,3,5...=1,3,5... a B 
16 1 
where = (42) 
™? mn 


The equations for the deflected surface, shear forces and bending moments can 
also be represented by double trigonometrical series. Each term of these series is of 
the form given in equations (22), (24), (25). 


The expressions for the bending moments only will be worked out in full. From 
(24) and (25) these can be expressed as 


oo 1 m2 b m X ny 
M,=- = —P { Wan (= re Yun )} 
m=1,3... Bt k B Do, k Do, 
b m Y 
M,=- { W (1.5 + (m3 Xn -n 3") x 


Substituting the values of Win, Ximns Ym» from equations (28), (29), (30) into (43) 
and (44) gives 


1 x m4 
+ (2 — + n* sin sin ‘ (45) 
+ n‘ Z.. sin —— sin = (46) 
where 
2 2 
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BUCKLING AND BENDING OF ORTHOTROPIC SANDWICH PANELS 


A corrugated core panel, with the corrugations running parallel to Ox, will now be 
considered. Substituting equations (34) and (35) and 


in equations (45), (46) and (47) gives 


1 x m? ir ( k 
>? k— +sr(1-: 
] x x m? u? u 
mn? | MEX. 


B? (i ata) tha + B? k(1+4) 


(51) 


The maximum values of the bending moments, M, and M, can be expressed as 
(M,) mex = 2,95" . ‘ (52) 
(My) max = . (53) 


where z, and z, are functions of - N,/P,, r, k, and 8. It can easily be shown that the 
position of maximum bending moments in the y-direction is on the centre line, y= 6/2, 
but the position in the x-direction has to be found graphically. Fig. 10 gives the bending 
moment (M,,),-»,2 along the centre line of the panel for the case of B=3-0, r=0-333, 
and k=1-25. For this case it is seen that the peak bending moment occurs at approxi- 
mately one-sixth of the length of the panel from either end. The values of z, and 2, 
for the case B=3-0 and k=1-25 are plotted in Figs. 11 and 12. The effect on z, of 
changing the aspect ratio @ of the panel is shown in Figs. 13 and 14. 


For a sandwich panel having zero core flexural stiffnesses, the bending moments 
are obtained by substituting K=1 and equations (34) and (35) in equations (45), (46), 
and (47):— 


1 x x mn? 1 


+ ont \ gin ; (54) 
a b 


mn 
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n* _.. MEX 


(56) 


+n? 


Hence the maximum bending moments of the panel can again be expressed in the 
form 


(Mz) max a,qb* (57) 


where a, and a, are functions of - N,/P,, r, ¢, and 8. 
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The Lift of Twisted and Cambered Wings in 
Supersonic Flow’ 


G. N. LANCE, MSc., Ph.D. 


(Department of Engineering, University of California) 


Summary: A generalised conical flow theory is used to deduce an integral 
equation relating the velocity potential on a delta wing (with subsonic leading 
edges) to the given downwash distribution over the wing. The complete 
solution of this integral equation is derived. This complete solution is 
composed of two parts, one being symmetric and the other anti-symmetric 
with respect to the spanwise co-ordinate; each part represents a velocity 
potential. For example, if y is the spanwise co-ordinate and x is measured 
in the free stream direction, then a downwash of the form w= - 2,, Ux|y| 
is symmetric and will give rise to a symmetric potential, whereas 
w= —a,, Ux|y|sgny is anti-symmetric and gives rise to an anti-symmetric 
potential. The velocity potentials of such flows are given in the form of 
Tables for all downwashes up to and including homogenous cubics in the 
spanwise and streamwise co-ordinates. Table III gives similar formulae in 
the limiting case when the leading edges become transonic; these are 
compared with results given elsewhere“) and serve as a check on = results 
of Tables I and II. 


1. Introduction 


In a recent paper“) the author described a method of obtaining the lift on a 
twisted delta wing, with subsonic leading edges. The method was applied to obtain 
the lift distributions in the cases when the downwash w over the wing plan form 
was proportional to x’y’, with O<r+s<3. However, the results obtained therein 
were confined to the cases for which the velocity potential, and hence the lift 
distribution, had the same symmetry, with respect to the spanwise co-ordinate y, as 
had the downwash. For example, if w= — ,, y then the velocity potential obtained 
was an odd function of y; if w= —«,,xy’ the lift distribution was an even function 
of y. The present paper gives the extension required to obtain the potentials when 
the symmetry of the downwash distribution is carefully displayed. 


For a delta wing with supersonic leading edges the velocity potentials have 
been calculated, for w= —,,x”y’ by Hunn and the author”. The symmetric and 
the anti-symmetric velocity potentials were obtained for each downwash distribution. 
Such distributions are useful in flutter and aeroelastic problems when it is desired 
to describe the wing shape, approximately, by means of a power series expansion. 
The expansion must be even in y for symmetric motion and odd in y for anti- 
symmetric motion. 


Received June 1954. 

*This work was carried out under the auspices of the F.O.A. Scientific Research Project 
TA 01-101-3006 (OEEC 151). 

National Academy of Sciences Research Fellow. 
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The method described in Ref. 4 is not ideal when both the symmetric and 
anti-symmetric velocity potentials are required, so first a modification of this method 
is described. The present method has been partially described by Lomax and 
Heaslet’ and by Lomax, Heaslet and Fuller’. To make the extension clear, the 
work of Refs. 6 and 7 is given, in outline, here. The lift distribution was obtained, 
in these two papers, as the solution of a certain integral equation. However, the 
solution given was not complete, a fact which explains why, in obtaining certain 
unknown constants, more equations were available than were required. (The same 
error was committed in Ref. 4). It is shown in this paper that there are, in fact, 
exactly the same number of equations as unknowns when the complete solution of 
the integral equation is used. (This point is discussed more fully in Section 3). 


NOTATION 
a= Btany 
ai, b,,c,,b,’,c,’ constants defined by equations (14), (15) and (20) 
Ap _ difference in pressure between lower and upper surfaces of wing 
t variable of integration 
u,v,w components of perturbation velocity 
x,y,Z rectangular Cartesian co-ordinates (Fig. 1) 
x,y, variables of integration 
B = (M?-1)! 
D,,...,D, functions of a defined in Tables I and II 
H(@,») kernel function, see equations (11) 


K,E complete elliptic integrals of the first and second kinds of 
modulus (1 — 


Mach number 
U _ free stream velocity 


@,, constants of proportionality for w given over the wing plan 
form, see equations (2) 


y  semi-angle of the delta wing 

n = By/x 

n, = By,/x 

«x degree of the downwash distribution 
p density of free stream 


6 variable of integration 


sgny = y/|y| 
BF [ denotes the “finite part” of the integral 
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LIFT OF TWISTED AND CAMBERED WINGS 


Fig. 1. 
The co-ordinate system and notation relevant to the delta wing. 


2. The Integral Equation 


The co-ordinate system is shown in Fig. 1. The assumptions made here are 
the same as those of Ref. 4. Attention is confined to an unyawed delta wing placed 
in a supersonic flow whose free stream velocity is U; linearised supersonic theory 
is used and the flow is of generalised conical type. 


The equations of the leading edges are 
which is written n= +4, 
where the non-dimensional co-ordinate 
n= By/x, 
a= Btany 
and B=(M?-1)}. 
Since subsonic edges only are considered, a< 1. 
The downwash, which is known over the wing, will be of the form 
w= —2,,Ux" | y|* (2a) 
for symmetric deformation of the wings, or 
w= —2,,Ux"|y|’sgny. (2b) 
for anti-symmetric deformations. 
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The local load coefficient can be written, according to linearised theory, 


Sp _ _ 4% (3) 


where p, and p, are the pressures on the lower and upper surfaces, respectively. 


Consider an elementary ray of the wing (Fig. 2) emanating from the apex O 
and inclined at an angle 6 to the x-axis. Suppose that the load carried by this 
element is 


where C is a constant. Then the upwash field due to a complete triangle with one 
side along y=0 and the other along y=x tan 83 is given by‘? 


(x—x,) Ap (x,, y,) dx, 


U 


where 7 is the area of the wing plan form that lies inside the fore Mach cone from 
the point (x, y,0). The generalised principal part of the integral is to be evaluated, 


as indicated by the symbol F 


Writing »,=By,/x and performing the integration with respect to x,, after 
using (4), we obtain 


(1+n) 
(= y U 4n | a n) dn, ; 0<a, l<n<a 
(6a) 
a(1—») 
=e" n,* y ya 
4n a "> n 
0 
(6b) 


with slight modification if a is negative. 


The element of upwash produced by the ray of angle dé is obtained by 
subtracting two such expressions. Let the upwash of the plan form with vertex 
angle 5 be 


w=f(a,n). . ; ‘ (7) 
of 
Then dw =f (a+ da, n)—f (a, 
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: 


x 
Fig. 2. 


The elementary ray, of apex angle d5, emanating from the apex of the wing. 


If this differentiation is applied to (6) and », is replaced by 


_ 
(a-t) 


there results 


(2) dw _ —BCda(x+1)a* dt 
4n t(a—t)*? 


4] dw _ De 
4x t(a—t)*? 


In one of these integrals the Cauchy principal part is required. 


» . 


-, -l<a<n<l.. 


(8) 


(9a) 


(9b) 


Hence the total downwash contribution from all such triangular elements lying 


between —a and +a is 


(2) 
where H (@,n)= ae 
at 
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By expressions (2a) and (2b) it is seen that the boundary condition to be 
satisfied by (10) is that Bxw/(x*U) is a polynomial of, at most, degree « in the 
variable || for symmetric deformations and a polynomial in || multiplied by 
sgn »(=m/|n|) for anti-symmetric deformations. Hence 


B(e+1) 


Performing the differentiation with respect to » gives 


(1 — 1?) (6) dé 


Thus the function C (6) is found to be a solution of the integral equation 


FA =3 for symmetric deformations . (14) 
=% a; for auti-symmetric deformations. 
i=0 


3. The Solution of the Integral Equation 


It is shown in the Appendix that the complete solution of (14) can be written 
in the form 


+ cosh-! 


and it remains to show how the constants b, and c, may be obtained. 


Since (16) 
(10) can be written 
1) (o 6*C (6) d0 | dt . (17) 


with a similar expression if » <0. 
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The order of integration in (17) may be reversed and a careful analysis 
shows that 


(1-1) 


t “> (f-a’)! 
(18) 


Now B«w/(x'U) is simply a power of || or a power of |»| multiplied by »/| 1]; 
thus, by equating like powers of this co-ordinate on each side of (18) it is seen that 
there are exactly the same number of equations and unknowns, once it has been 
decided whether the odd or even solutions are being sought. For example, if «=2, 
there are just three equations and the six unknown constants b,, b,, b., b3, c,, Co. 
In seeking, say, the anti-symmetric load distribution only the constants with odd 
subscripts are selected, and the remainder are equated to zero. Thus the unknowns 
are reduced to exactly three. 


4, Examples 


Two examples are given which bring out the technique for obtaining the 
constants. As a final result the velocity potential ¢,, is required for downwash 
distributions given by equations (2). Since »=By/x, dx=—Bydn/n?, and thus, 
by (3), 


Ap __ 4 % 


It follows that we must assume an expression for ¢ of the form 


k+l a + 
g=aU xt} E cosh- +3 b,’ 1’ ‘ (20) 
1 0 
which, when substituted in (19), will give an expression for 


&Ap 
$pU*y* 


which is of the correct type, namely expression (15). 


(A= (21) 


For illustrative purposes, only the cases 
(a) w= —a@,,Ux|y| 
and (b) w=-—a,,Ux|y|sgny, 


are considered, in which case kx=2. Case (a) corresponds to the symmetric and 
case (b) to the anti-symmetric velocity potentials. 
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Case (a) 


From (20), the expression for ¢ is 


= +b,’ +C,'n? cosh-! | (22) 


where only the coefficients of even powers of 7 have been selected and the factor 
B-? is inserted to simplify b,’, b,’ and c,’. 


Now (19) and (21) show that 


(23) 


Comparing this equation with (15) shows that 


12a’b,’, 
b,=4(a’b,’ +ac,’-2b,) 
and 
Thus (18) becomes 


t t 


442 -1 
t? cos 


1 


By equating the coefficients of 1 on each side, the following three equations are 
obtained : — 


(3a*b,’ +t? {a*b,’ + ac,’—2b,'}) 


(t? — a’) 


1 


a 
] dt=1 (25) 


F 442 -1 
fa t?) cos 
1 


1 
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With the substitution ¢? =a? nd’ u equations (24) and (26) reduce to 


ab,’ [2a*K (1 + a?) E)+ ab,’ [a’ K (3 
—c,’ [2a?K (3—5a?+a*)— (27) 


ab,’ (5 ab,’ [2a*K—(a’ +a*) E]+ 
+a’c,’ [a’K (1+ a’)—2a@’°E]=0 (28) 


respectively, while (25) can be evaluated in terms of elementary functions and gives 


The tables of integrals of elliptic functions by Byrd and Friedman were found 
most useful in all these evaluations; throughout this work K and E denote the 
complete elliptic integrals of the first and second kinds, respectively, whose modulus 
is (1—a*)}. 

After substituting for c,’ from (29) in (27) and (28), we solve for b,’ and b,’ 
and obtain 


b,’=a[3a' K*— KE +(4+ a’) E’]/D, G0) 
and b,’= —2a[a*t K? (1+a’) KE+3a°E’]/D,, & 
where D,=Sa' K? 8a? (1 +a’) KE—(4— 19a?+ 4a‘) E’; 


this completes the derivation of ¢,, of Table I. 


Case (b) 


The anti-symmetric solution of (20) is now required, namely 


From (19) and (21), we have 


b,’ (2a’ _1|4@ a 


ai 
(34) 
and (18) gives, in this case, 
c,/at’ 
(35) 


Equating coefficients of 7 on both sides of this equation, we have 
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t (t? — (f-a) 
(36) 
1 
(37) 
f war a a b,’t (2a? 
(38) 


« In the present case (36) and (38) can be evaluated in terms of elementary 
functions and they show that 


With the elliptic function substitution t?=a’ nd? u, (37) gives 


' [a +a") 


This completes the derivation of ¢,, of Table II. 


(40) 


5. Results 


The results for the symmetric velocity potentials are listed in Table I. They 
are valid for -a<n<a. The velocity potentials which are anti-symmetric in 1 
are listed in Table II; these also are valid for -a<»<a. (In all tables || is used 
in the argument of the cosh~' terms). It is apparent from these complicated results 
that their derivation is long and tedious. A check is desirable; it is applied by 
taking the limit as a tends to 1 from below, and comparing these much simpler 
formulae with the corresponding limiting values derived from the supersonic 
formulae“? by letting a tend to 1 from above. In every case the results agree, so it is 
unlikely that any errors are now present. 


6. Conclusions 


The formulae listed in Tables I and II give the symmetric and anti-symmetric 
velocity potentials #,,, respectively, when the downwash over the wing plan form is 
given in the form w=-—2,,Ux'y’. Table III gives the formulae for wings with 
transonic leading edges. 


The Tables can be used for flutter calculations as described in Ref. 5, Table | 
being used for symmetric modes of vibration and Table II for anti-symmetric modes. 
The formulae given herein can also be used in aeroelastic calculations, the shape 
being specified in a power series of the form =2,,x"y’. The technique used is 
similar to that of Ref. 3 for swept wings with moderate aspect ratio. 
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The method presented can, in principle, be extended to values of « > 3, but the 
algebraic work gets more and more laborious, and it is doubtful whether any 
useful purpose would be served by carrying this work further. A slight extension 
only would be required to include thickness effects, and this will form*the subject 
of another paper. 
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Appendix 
It is required to show that 
la| 
& cosh-? | —|. 
f(@) t 2% cos (41) 
is a solution of the integral equation 
= a; | for symmetric deformations. (42) 
i=0 


= a,|n|' for anti-symmetric deformations. 


The first summation in (41) is the usual solution and the fact that it satisfies (42) 
will not be proved. 
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Consider the expression 


cosh" || a9 a cosh~" dé « cosh-? d@ 


6 
1 2an cosh-* di 
0 
1 1 
nia cosh~? dt 
0 nla 
Integration by parts gives 
nla d 1 
t t dt 
0 n/a 
1 1 
t dt 
{. n+a: t(1 t(1 — 
0 n/a 
These two integrals can be found in tables (e.g. Ref. 1) and we get 
{ sin =-4r?—, if » 0. . (44) 


If »=0, (43) shows that J,=0. The exact value of J, is unimportant in this context, the 
important point being that it is a function of »/| | only. 


k=0 
_1|4 
6" cosh \5 


-a 


the expression in { } is clearly a polynomial of degree r in || or a polynomial of 
degree r in || multiplied by »/|»| according as the motion is symmetric or anti- 
symmetric. Hence it follows that the second summation in (41) gives a contribution, 


to the integral equation (42), of the form for symmetric, or | n|‘(y/| 1!) 
1 1 
for anti-symmetric motions, which is thus of the correct form. 


Thus it has been shown that the complete solution of the integral equations (13) 
or (14) may be written in the form (15), despite the fact that the left hand side of (10) 
and its derivatives may be discontinuous at 7=0. It has been shown, by Lomax’, that 
the solution of (13), having a discontinuity at 9=/, where -a<1<a may be derived 
from the standard Séhngen inversion formula. 


*Private communication to the author. 
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TABLE I 


EXPRESSIONS FOR THE SYMMETRIC POTENTIALS 
$,s i8 the velocity potential for the flow whose downwash is of the form 


w= —a,,Ux"|y |*. 
Ux (a? 
Ux?a (E—a?K) (a? — | 
a, Ux? (1 —a?) (a? — 


2a, 


Ux*{a? —(1 + a2) E]— 7? (5 —3a?) + 2E (1 — Sa? + 3a*)]} (a? 


2 


3B°D, 


{(3a*K? — —8a°KE + (4 +a?) E?—2n? [a*K?— 2a°KE (I +a?) (a? — 
7B? D, 


1, Ux*n? cosh-? 


| 
+ 


7B? 
Ux' {a°K (3 — 5a?) —2E (3 —5a? +a*)+ n? [2a7K —E (1 (a? — 


3BD 


1 


where D, =5Sa*K? —8a°KE (1 + a”)—E? (4— 19a? + 4a*) 


a? [2a*K? (3a? —7) + 4a?KE (8 — a? — 3a*)—2E? (4+ lla? —11a‘)]+ 


2 2 4 6 | 2 
latK? (27—55a? + 48at)—a*KE (48 — 58a? — 46a‘ + 96a*) — (a? — + 
| —E? (12—119a? + 239a4 — 152a®)] 
Uxtn* 
| 
27B4 
2 a,, Uxt (1—a*) [a?K (9 —a?)—E (3 + 7a? —2a*)] — } (a? — 12)! 
Pre n? (27 —31a? + 12a*)+E (6 — 55a? + 65a* —24a‘)] 
2 [a*K? (78 — 88a? — 13a* + 3a*)—2a2KE (96 — 98a? — 16a* — 5a* + + 
Uxta | +E? (72— 10a? —84a*—9a® + 11a8)] + 9? [a*K? (306 — 51 1a? +41a* + 
| + 48a*)—2a2KE(312 + 10a? — 559a* + 73a°+48a*)+ E? (244+974a? — 7 
— 1107a*— 159a* + 152a8)] 
a,, Ux*n? [12+(2 —a?) n?] cosh-? 4 
2a,, Ux* (1—a*) 
+1? [a?K (9—a”)—E (3 + 7a? —2a*)] 7 
where 


May 1955 


D, =a*K? (27—31a? —8a*) — 2a*KE (24 — Sa? — 15a* — 16a*) — 
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TABLE II 
EXPRESSIONS FOR THE ANTI-SYMMETRIC POTENTIALS 


#, is the velocity potential for the flow whose downwash is of the form 


w=—a,,Ux" | y |* sgny. 


Uxn cosh-! | 


B? [a?K —(2 —a?) E] 


2a,, Ux? cosh~! 
—2a,, Ux?an [a*K —E] (a? 


[a2K —(2—a?) E] 7B 


2 Ux'n® cosh-* F | 
2 a. Ux%an (1 —3a?)—2(1 — 2a) E] (a? — 


3xB° 
a, , Ux? (1 —a?)?n (a? — 
B?D 


3 


a4, (6 +n?) cosh-? F | 
ay, Ux*an [a?K (23—9a%)— 2E (8—a2)] 


3nBD, 


where D, (1+a?)—2E (1 —a? +a‘). 


a*[a?K (1 —9a?)—E (2 —7a? —3a‘*)] + 


4 1— 2 


BD, 
Uxtan (8-570? + 45a*)— 2a2KE (16 —89a? + 69a4) + 
+E? (32 — 128a? + 83a* + 9a*)+2a?n? [a*K? (7 —3a?)— (a? — + 
—2a°KE (8 — a? —3a*)+ E? (4+ 1a?—11a*)] 
2a,, Ux*y? cosh- | 
" 
(12—31a? + 27a*)—(24—65a? + 5Sa* —6a*) E + 
4(1—@?2 
3B°D, 
a4, Uxtay { (160 18702 + 63a*)— 2a*KE (200-2230 
B2 (160 —84a? — 31a*—9a®) +2 (1 —a?) n? [a*K? (13 + 3a) — | + 
| —2a*KE (16 + 5a? +3at) +E? (12-+13a2+7a*)] | 
a 
Ux*n (2 +n?) cosh? 
+ 
7B 
where 


D ,=a*K? (8 + 31a? —27a*)—2a°KE (16 + 15a? + Sa* — 24a°) +B? (32 — 64a? — 151a* — 119a° + 
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LIFT OF TWISTED AND CAMBERED WINGS 


TABLE Ill 
EXPRESSIONS FOR THE POTENTIALS WHEN THE LEADING EDGES ARE TRANSONIC 


(b) Anti-symmetric 


4a,, Ux?n 


Uxn cosh-! 


2a,, Ux?» cosh~! 
2&9 


3nB 


2 (1-7) 


1 
2 Ux*n* cosh-! | 


16 a, , (l1—n?} 
157B* 


(a) Symmetric 
Ux 
Br 
1 
& 4, Ux?n? cosh-? 
& 4, Ux? 
Por 3nB? TB? 
4a,, Ux? (1—n?) 
Pro 3nB 
Pos 457B° 
a, , Ux? (3—2n?) (1—n?) 
15xB? 1 
a,, cosh-?| 
+ 7B? 
4a,, Ux* (11+?) (1—n?) 
P20 
707 B*4 1 
Pos a, Uxtn* | — | 
27B4 
16 Uxt (1+17n?) (1 —7*) 
Pi 
Ux* (54—125n?) (1 — 
4207 B? 
1 
(12+?) n? cosh=? -| 
+ 
127B? 
6 16 a,, Uxt (5+?) 
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Ux*n (6+?) cosh-? | 


+ 


1 


3nB 


16 Bos (i+ 5?) (1 
1057B* 


2 Ux*n (15 —2n?) (1 —n?)# 
2a,. Ux*n’® cosh-! 


Ux*n (17 +n?) (1—1n?)+ 
3157B? 


105nB 
(2+ n*) cosh-! 
+ 


= 
| 
+ 12a°*) 
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REVIEW 


Review 


Studies in Mathematics and Mechanics, presented to Richard von Mises by Friends, 
Colleagues and Pupils, and edited by Garrett Birkhoff, Gustav Kuerti and Gabor 
Szegd. New York Academic Press, 1954. 353 pp. $9.00. 


This tribute to the eminent work of Richard von Mises in applied mathematics and 
mechanics takes the admirable form of a series of 42 articles, on branches of these 
subjects in which he had worked, by “ friends, colleagues and pupils.” It is preceded 
by a valuable bibliography of von Mises’s writings, and a brief but penetrating study by 
Philipp Frank of the views and aims which governed the selection by von Mises of a 
wide variety of scientific, mathematical, philosophical and literary subject-matter for 
investigation. Here are noticed briefly only four articles in the volume, which seem 
to be those of most interest to aeronautical scientists. 


Professor Temple contributes an article on his very important concept of “ weak 
«functions ”” or “ generalised functions,” which is treated further in J. London Math. 
Soc. 28, 134 (1953) and Proc. Roy. Soc. 228, 175 (1955). In these investigations he 
gives a simple and rigorous foundation to the theory of delta functions, principal values 
and finite parts of integrals, and the Fourier analysis of singular functions of all kinds. 
Much of this work greatly simplifies the analysis in many well-known kinds of aero- 
dynamical calculation. 


A simple model of shear flow turbulence based on considering the motion from a 
Lagrangian point of view and postulating certain resistances and rotations for fluid 
elements depending on their difference of velocity from the mean and on the local mean 
shear is put forward by Dr. J. M. Burgers. 


The oscillating vortex wake behind an obstacle at Reynolds numbers of order 10? 
is studied by means of the Oseen approximation in a paper by Professor C. C. Lin. 


Professor Howard Emmons derives a general correlation for experimental results 
on natural convection, film condensation, film boiling and film melting, by re-analysing 
the significance of the relevant Rayleigh and Reynolds numbers. 

M. J. LIGHTHILL. 


Errata 


It is regretted that in the paper by L. E. Fraenkel “ On the Unsteady Motion of a 
Slender Body through a Compressible Fluid” published in the February 1955 issue of 
The Aeronautical Quarterly (Vol. VI, p. 59), the printed impression of some of the pages 
was very poor. 


The most serious result of this was that the identity between equations (16) and 
(17), page 65, which should read 


apparently contained “ z” and not “ z?” in the expression on the right hand side. 


In addition the letter “ d” should have appeared before the “ 6” at the very end of 
equations (76) and (77) on p. 80. 
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